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Abstract 



How much information about the shape of an object can be inferred from its image? 
In particular, can the shape of an object be reconstructed by measuring the light it 
reflects from points on its surface? These questions were raised by Horn [HO70] who 
formulated a set of conditions such that the image formation can be described in terms 
of a first order partial differential equation, the image irra.dia.nce equation. In general, 
an image irradiance equation has infinitely many solutions. Thus constraints necessary 
to find a unique solution need to be identified. 

First we study the continuous image irradiance equation. It is demonstrated when 
and how the knowledge of the position of edges on a surface can be used to reconstruct 
the surface. Furthermore we show how much about the shape of a surface can be 
deduced from so called singular points. At these points the surface orientation is 
uniquely determined by the measured brightness. 

Then we investigate images in which certain types of silhouettes, which we call 
b-silhouettes, can be detected. In particular we answer the following question in the 
affirmative: Is there a set of constraints which assure that if an image irradiance 
equation has a solution, it is unique? To this end we postulate three constraints upon 
the image irradiance equation and prove that they are sufficient to uniquely reconstruct 
the surface from its image. Furthermore it is shown that any two of these constraints 
are insufficient to assure a unique solution to an image irradiance equation. Examples 
are given which illustrate the different issues. 

Finally, an overview of known numerical methods for computing solutions to an 
image irradiance equation are presented. 
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Chapter I 



Motivation 



How much information about the shape of an object can be inferred from its image? 
We are interested in a special aspect of this question: the reconstruction problem, which 
is to determine the shape of an object from measurements of the light reflected from its 
surface. In general, there are many surfaces which can give rise to the same image. So, 
we will try to identify and analyze constraints such that the shape of a surface can be 
uniquely reconstructed from its image. Our search for these constraints is guided by the 
following question: Are there any "properties" of an object that are also "properties" 
of its image or vice versa? In fact, as we shall see, such properties do exist. 

Our work is based on Horn's thesis [HO70]. He formulated a set of conditions (to be 
discussed in the next chapter) which lead to a relation between the perceived brightness 
of a small patch of a surface and its normal vector. This relation, the image irradiance 
equation, is a first order partial differential equation (abbreviated in the following by 
FOPDE) and each of its solutions determines the shape of an object. The problem of 
finding solutions to the image irradiance equation is referred to in the literature as the 
shape from shading problem. 

We will take two approaches towards finding a solution to the shape from shading 
problem termed as the local and the global approach. By the local approach we mean 
that only a small patch of an image is used to determine the shape of a surface. To 
the contrary, in the global approach we examine images in which a silhouette can be 
detected (here we refer to the outline of an image as a silhouette). 

Intuitively, it seems clear that by looking at an image in which a silhouette can be 
identified we should be able to conclude more about the shape of a surface whose image 
we are analyzing than by just looking at a little patch. Wc will show that from certain 



images which contain a silhouette we can uniquely infer the shape of the surface which 
gives rise to that image. Unfortunately the global approach is not always satisfactory; 
there are also many images containing silhouettes which could be the images of infinitely 
many different surfaces. There are also infinitely many surfaces which locally look the 
same. So we will determine conditions under which the global approach is better than 
the local approach. Notwithstanding, one can sometimes draw interesting conclusions 
about the shape of surfaces which give rise to the same image by just looking at a small 
patch of this image. 

The local approach is taken to an extreme when we pose the following question: 
What can be deduced about the shape of a surface from so-called singular points of an 
image irradiance equation? At these points the surface normal to all solutions to such 
an equation is uniquely determined by the brightness there. We investigate the above 
stated question for a certain class of image irradiance equations, the so-called eikonal 
equations, which describe a variety of physical phenomena. For instance, experimental 
data suggest that the flux of secondary electrons in a scanning electron microscope can 
be described by an eikonal equation [LAWH60]. By using these secondary electrons to 
modulate the appropriate devices, an image of a surface is created by the microscope. 
Such an image exhibits shading [HO70, pp.85-87] and therefore to determine the 
shape of a surface from its image one effectively has to solve an eikonal equation. 
Studying eikonal equations, we show that the absolute value of the Gaussian curvature 
at a singular point of all surfaces which give rise to a particular image, is the same. 
Furthermore, assuming that the surface is convex at a singular point, we show that 
its shape can be uniquely determined in some neighborhood of such a point from the 
image intensities alone. 

The other aspect of the shape from shading problem which we explore is its solution 
when the image contains a b-silhouette (which is defined below). In this case a global 
approach is taken. Let us first define the bounding contour of a surface: a point P 
is on the bounding contour if the line connecting the viewer and P grazes the surface 
(i.e., if this line lies in the tangent plane of P). Furthermore we assume that no two 
parts of a surface obscure each other, i.e., we assume that the bounding contour is not 
an occluding contour. The image (assuming orthographic projection) of a bounding 
contour will be called the b-silhouette. The surface normal at a point on a bounding 
contour is parallel to the normal vector to the b-silhouette and both vectors lie in the 
same plane. Thus, some or all of the first order partial derivatives of the function 
defining the surface are infinite for points on the bounding contour (we will say that 
some components of the surface gradient are singular along a curve). 

Remark: In the context of this report the phrase a function f(x,y) is singular at a 
point (io, t/o) will always mean that: 

lim /(z,y) = ±oo. (I-l-l) 

X— +XQ 

v-»vo 

This terminology should not be confused with the notion of a singular solution which 
we will explain in section A.4. 
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. For example, the bounding contour of a hemisphere lying on a plane B is a circle. 
Consider a Lambertian surface which has the property that each surface patch appears 
equally bright from all viewing directions. If we look at a Lambertian hemisphere such 
that the viewer and the light source are at the same point, its b-silhouette can be 
determined from its image. In this case the image irradiance equation describing the 
imaging situation is singular and all the surfaces which satisfy such an equation have 
a bounding contour. (We extend the notion of a singular function to equations.) 

On the other hand, the existence and the position of a b-silhouette cannot be 
determined from a continuous image irradiance equation. Thus the existence of a 
surface which satisfies a continuous image irradiance equation and which has a bounding 
contour does not necessarily imply that all surfaces which satisfy this equation have a 
bounding contour. 

In this report, Horn's work on the shape from shading problem is extended. 
He studied primarily images of smooth surfaces where the imaging situation can be 
described by a continuous image irradiance equation (which is defined rigorously in 
section III.l). Yet some objects have edges and their surface normals assume different 
values depending upon which side an edge is approached from (we will say that the 
surfaces normals are discontinuous along an edge). Can an edge be invisible in an 
image? In other words, is it possible that a surface has an edge without it producing a 
discontinuity in the equation? As it turns out, discontinuous solutions can arise from 
continuous equations. In this case initial conditions (discussed in sections III.2, A.6 and 
A.7) provide information about the occurrence of a discontinuity. Additionally, we will 
show (in section III.2.3) under which conditions edges on a surface can be used as such 

initial data. 

We keep in mind that the image irradiance equation describes a physical situa- 
tion and therefore the only solutions considered are those corresponding to (piecewise) 
smooth surfaces. (A rigorous definition of the notion smooth surface is given in section 
III.l.) However there are image irradiance equationslor which no solution correspond- 
ing to a smooth surface exists. In particular, the existence of such a solution to a 
singular image irradiance equation is not guaranteed. We are primarily concerned with 
the identification of constraints which allow one to solve the reconstruction problem 

uniquely. 

In general, a PDE describes a class of processes rather than a particular one. 

Consider, for example, the Laplace equation: 

A/ = (1-1.2) 

where A denotes the Laplace operator and / a scalar field. This PDE constrains the 
sources and the curl of the field / to be zero, but an infinite number of different fields 
exist which fall into this category. Only when some further conditions about / are 
specified can the solution to the Laplace equation be restricted to a single one. 

Similarly, there are an infinite number of different surfaces which satisfy a given 
image irradiance equation. Thus, as Horn [HO70, II075] has already observed, in 



general the image irradiance equation alone is not sufficient to solve the reconstruction 
problem uniquely. It remained an open question whether there are any imaging situa- 
tions for which every surface gives rise to a different image. As indicated before, taking 
the global approach towards finding a solution to the shape from shading problem 
allows us to answer this question affirmatively. To this end we will postulate three 
constraints upon the image irradiance equation and prove that if these constraints are 
known to hold, the reconstruction problem can be solved uniquely. 

We now briefly describe the different chapters in this thesis. 

Chapter II is a summary of the issues involved in the shape from shading problem. 
For a more extensive discussion of this material see [HO70], [H075] and [WOOD78]. 

Chapter III gives an overview of the different mathematical problems involved in 
solving a FOPDE. The mathematical details used in this chapter can be found in 
appendix I. In section III.l we define two classes of image irradiance equations, the 
continuous and the singular equations. Section III.2 deals in detail with the continuous 
image irradiance equation. In particular we exhibit in section III.2.2 that certain 
continuous image irradiance equations can be transformed into singular equations. 
This implies that all surfaces which satisfy such an image irradiance equations have a 
bounding contour. In section III.2.3 we examine how edges on a surface can be used 
to reconstruct the surface. To help us understand the variety of integral surfaces of an 
image irradiance equation we introduce in section III.2.4 some concepts from dl ff eren ™ 
geometry and in section III.2.5 gradient space as popularized by Mackworth [MAC73J 
and Horn [H077]. Section III.3 describes the basic issues involved in solving a singular 
image irradiance equation. In section III.3.1 Marr's [MA77] work on occluding contours 
is reviewed. Section III.3.2 discusses the method of characteristic curves for singular 
image irradiance equations. 

In chapter IV we show how a singular point of an eikonal equation constrains its 

possible solutions. 

In chapter V we postulate three constraints upon an image irradiance equation sucn 
that the reconstruction problem can be solved uniquely. In section V.l it is proven that 
if an image irradiance equation satisfies these constraints it has a unique solution. It 
is shown in section V.2 that any two of these constraints are insufficient to assure a 
unique solution to an image irradiance equation. 

Chapter VI gives a brief description and discussion of known numerical methods 
for computing solutions to an image irradiance equation. 

Chapter VII summarizes the results of this report and suggests some possible 

applications. 
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Chapter II 



The Shape from Shading Problem 



There are basically three components to the shape from shading problem which 
have to be taken into account. They are the light source, the object and the camera as 
depicted in figure 1 which is taken from [WOOD78, p.32], and termed as an imaging 
configuration. Henceforth we will assume that an image of a surface is produced by 
a camera. The shading of such an image can be explained as follows: The exposure 
of film in a camera (for fixed shutter speed) is proportional to image irradiance, the 
light flux per unit area falling on the image plane. Similarly, grey levels measured in 
an electronic imaging device are quantized measurements of image irradiance. It can 
be shown that image irradiance m turn is proportional to scene radiance, the light flux 
emitted by the object per unit projected surface area per unit solid angle [HOS79]. 
The factor of proportionality depends on details of the optical system, including the 
effective /-number. 

Scene radiance depends on the 



• surface material and its microstructure, 

• the incident light flux, and, 

• the orientation of the surface. 



Now we want to relate the shape of a surface to the shading of its image. Consider 
a viewer- oriented coordinate system with the viewer located far above the surface on 
the z-axis. If the objects imaged are small compared to their distance from the viewer, 
one can approximate the imaging situation by an orthographic projection: 
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SURFACE 
NORMAL 




Figure 1. Imaging configuration 
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Zq 



(II.l.l) 



where (z,y) are the coordinates of the image of a point {x,y,z) made with a system 
of effective focal length /, and the viewer is at a distance z above the origin. We 
assume that (x 2 + y 2 + z 2 ) < z\. For simplicity and without loss of generality it is 
also assumed that the viewing direction coincides with the z-axis. 

The orientation of a patch of a surface can be specified by its gradient (p, q, —1), 
where p and q are the first order partial derivatives of z with respect to x and y. For a 
given surface material and known incident light flux, scene radiance will depend only on 
surface gradient. The function which describes this dependence, R{p, q), (or a contour 
representation in gradient space), is called the reflectance map. 

Recall that image irradiance and scene radiance are proportional and that we 
assume orthographic projection. If E{x, y) is the observed image irradiance at the point 

(x, y) in the image then: 

K y) R( P ,q) = E(x,y) (H.1.2) 

where (p, q) are two components of the gradient at the corresponding point on the object 
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being imaged. This equation is called the image irra.dia.nce equation. It is clearly a first 
order partial differential equation since it involves only the first order partial derivatives 
p and q and the coordinates* i and y. In summary, to derive an image irradiance 
equation we have to know the reflectivity function and the geometry relating the light 
source, the object and the camera. 

How can an image irradiance equation be used to analyze images? Informally, sup- 
pose enough information about the imaging situation is known so that the reflectance 
map can be derived. Then at every point (x ,yo) the image irradiance denoted by Eq 
can be measured. Using the image irradiance equation we can determine the possible 
values for p and q such that: 

R(p,q) = E . (H.1.3) 

However, the values for p and q cannot be chosen independently as p and q have to be 
the first order partial derivatives of a function z = z(x, y) defining a smooth surface. In 
general there are many values for p and q which satisfy equation (II. 1.3) and represent 
the components of a surface gradient. Consequently, an image irradiance equation has 
many solutions or equivalently, for a fixed imaging configuration many surfaces will 
give rise to the same image. Thus we shall determine constraints under which the 
solutions to an image irradiance equation are restricted to a unique one. 

A word of caution: Several issues such as mutual illumination, shadows or specularity 
are not addressed here. 
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Chapter III 



Overview 



111.1. Classification of the Image Irradiance Equation 

In this chapter, we will further investigate the shape from shading problem. Our 
interests are twofold: 

1) We wish to find the solutions to an image irradiance equation. 

2) We wish to determine constraints which guarantee a unique solution to an 
image irradiance equation. 

The results concerning these two issues are different for continuous and singular image 
irradiance equations, both of which are defined later in this section. First we make 
some general observations concerning an such equations and their solutions. 

An image irradiance equation is a first order partial differential equation in two 
variables x and y. We are looking for a solution z = z(x, y) (also called integral surface) 
which is a function of x and y and defines a surface from which light is reflected. 
To be precise, the class of functions which we allow as solutions to the FOPDE has 
to be specified and we proceed now with some relevant definitions. Using standard 
nomenclature, a function f(x, y) is said to be of class C k if it has continuous A;-th order 
partial derivatives. The following definition captures formally our geometrical intuition 
about a surface in ?R 3 : 
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Definition: Let B C 5R 2 be a connected and closed region, and let z.B *-+ 3fc 3 be 
continuous on B. Then the single- valued function z — z(x, y) defines a smooth 
surface if z is C 1 at every interior point of B. 

For smooth surfaces the unit surface normal is well defined for every point on the 
surface by: 

(i.o,ft)x(o,i,ft). m x 1} 

ll(i,o,ft)x(o,i,fc)|| 

A piecewise smooth surface is defined similarly. Unless otherwise stated we will always 
assume that the solutions to any given FOPDE are (piecewise) smooth surfaces. We 
will denote by p and q the first order partial derivatives of z = z{x, y) with respect to 
x and y. 

We shall exploit the following two facts about an image irradiance equation: 

1) The equation does not depend explicitly on z. 

2) The equation involves only two functions R and E, such that R depends 
only on p and q and E depends only on x and y. 

An immediate consequence of 1 is that there are no singular solutions to the equation, 
i.e., the complete integral describes all possible solutions. (For an explanation of these 
terms see section A.4). As for 2, we will study image irradiance equations which fall 
into either one of the following two categories: 

* 

1) The functions R and E are C 1 . We say that such an image irradiance 
equation is continuous. This case is discussed in section Hl.2. 

2) The function RisC 1 . The function E is singular in x and/or y but for all 
points (x,y) at which E(x,y) assumes finite value, it is C 1 . We say that such 
an image irradiance equation is singular. This case is discussed in section III.3. 

Since the measured image irradiance is always finite, it seems at first that singular 
image irradiance equations are not of practical interest. However, for our purposes we 
can think about such equations as useful mathematical constructs, i.e., a singular image 
irradiance equation can be obtained by transforming a continuous image irradiance 
equation appropriately. By appropriately we mean that the transformation is one-to- 
one and onto and that the solutions to the original equation and the transformed one 
are the same (section 1II.2.2). For example, there is a transformation between the 
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continuous image irradiance equation (III.1.2) and the singular equation (ffl.1.3): 

i + p 2 + <r 

D 2 , 2 = * 2 +y 2 , (in.i.3) 

We still have to explain why such a transformation is useful, i.e., why it makes 
more sense for us to analyze the singular equation (III.1.3) instead of the continuous 
equation (III.1.2). The reason is that we can gain some information about the integral 
surfaces of a singular image irradiance equation without first solving it. To this end let 
us define the terms bounding contour and b- silhouette: 

Definition: A point P lies on the bounding contour of a surface defined by 
z = z{x, y) if the tangent plane at P is perpendicular to the x-y plane. The 
b-silhouette is the orthographic projection of the bounding contour onto the 
x-y plane. 

We will show below that each integral surface of any singular image irradiance equation 
has a bounding contour, a fact which is not true for continuous image irradiance 
equations. In section III.2.2 we will construct two integral surfaces of a continuous 
image irradiance equation such that one of them has a bounding contour whereas the 
other does not. Why do all integral surfaces of a singular image irradiance equation 
have a bounding contour? Recall (chapter I) that for points on the bounding contour 
p and/or q become infinite. Now, an image irradiance equation is singular if there 
are points (z ,yo) such that the only values for p and q which satisfy the equation at 
(xo.Vo) are infinite. This explains the previously stated question. Furthermore, as the 
b-silhouette is the projection of the bounding contour onto the image plane, the points 
which constitute the b-silhouette can be uniquely determined from a singular image 
irradiance equation (section III.3). 

We will not be concerned with discontinuous reflectance maps since these are rare. 
They may occur when dealing with specularities, an issue not addressed in this report. 
(It is also assumed that the reflectance map is not a constant since then the image will 
also be constant.) 

III. 2. The Continuous Image Irradiance Equation 

In the mathematical literature, the most studied FOPDE's are continuous and 
have continuous first order partial derivatives. An overview of known results in this 
area can be found in appendix I and we summarize them in the next few sections 
only very briefly. In particular we explain why, in general, additional information is 
needed to restrict the solutions to a given image irradiance equation to a single one. 
Furthermore we address the problem of whether one can deduce any properties of the 
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integral surfaces of a given image irradiance equation by just inspecting the equation. 
To this end we investigate singular points, at which the tangent plane to all integral 
surfaces of an image irradiance equation is uniquely defined by the image intensity 
(section III.2.2). We can. then show (chapter IV) that for a class of image irradiance 
equations the absolute value of the curvature of the integral surfaces at a singular point 
is uniquely defined by the measured brightness there. In addition, we discuss a class 
of image irradiance equations for each member of which we can deduce that all of its 
integral surfaces have a bounding contour. 

The main tool for solving a FOPDE is the construction of so called characteristic 
curves. The following is a short explanation of this notion. A more detailed exposition 
can be found in sections A.2 and A.3. Let 

R(p,q) = E(x,y) (IH.2.1) 

be a continuous image irradiance equation. What kind of information about an integral 
surface (defined by z = z(x, y)) of the previous equation can be deduced from it? For a 
given point P on z the quantities p and q are constrained by (IH.2.1) to lie on a curve. 
Since p and q determine the normal (p,q, — 1) at P, equation (IH.2.1) constrains the 
feasible tangent planes at P to a one-parameter family, "enveloping a conical surface 
with P as vertex, called the Monge cone" [COHI62b, p.75]. The directions of the 
generators of a Monge cone are called characteristic directions. Now, the characteristic 
curves are those curves on an integral surface which at every point have as their tangent 
direction a characteristic direction. The characteristic curves can be determined by 
solving the characteristic equations which are five ordinary differential equations whose 
solutions depend on initial values: 

%--(pF. + F.) g—W. + F,). 

Equivalently, we can say that for distinct initial values, different characteristic curves 
are determined. Since each integral surface of a FOPDE is swept out by characteristic 
curves (which we prove in section A.3), initial conditions (which are discussed further in 
the next sections) are necessary to restrict the solutions to an image irradiance equation 
to a single one; in the case where an image irradiance equation is continuous, a unique 
solution to the problem can be found provided that an appropriate initial strip defined 
by x = .x(i),y = y{t),z — z(t),p = p(t) and q = q(t) is known (section A.6 and 
A.7). By appropriate we mean that this strip is not a characteristic strip and that the 
following condition holds: 

A ^dxdy_dydx^ } 

ds dt ds dt 
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Depending upon the initial curve, the integral surface is either continuous or discon- 
tinuous as will be discussed in greater detail later on. 

111.2.1. The Linear Continuous Image Irradiance Equation 

We first review the case where the image irradiance equation is a linear FOPDE 

e g ' 

ap + bq = E{x,y) (HI.2.4) 

where a and b are constants. This equation is of practical interest as its linear reflectance 
map describes the reflectivity properties of the maria of the moon where the constants 
o and b define the direction towards the light source (i.e., the sun) [HO70]. Linear 
FOPDE's are special cases of quasi-linear FOPDE's (in which a and 6 are functions of 
x, y and z). However, an image irradiance equation cannot be a quasi-linear FOPDE 
unless it is a linear FOPDE. Furthermore we can show the following lemma: 

Lemma: The solutions to an image irradiance equation of the form: 

f{ap + bq) = E{x,y) (HI.2.5) 

where / is a bijection, f- x {E{x,y)) is C 1 and where a and 6 are constants, 
can be obtained by a simple coordinate transformation from the solutions to 
an image irradiance equation of the form: 

p = E(x,y). (IIL2.6) 

Remark: A bijection is a one-to-one and onto function. In the above lemma, f~ 
denotes the inverse function of /. 

Proof: Since the function / is' a bijection, equation (III.2.5) and the transformed 

equation (III.2.7): , ttt 

ap + 6«= S /- 1 (JB(x,y))' (IH.2.7) 

have the same solutions. We abbreviate f- l {E(x,y)) by E{x,y). To prove the lemma 
we distinguish three cases depending upon the values of a and 6. 

Case 1) a ^ and b ^ 

The image irradiance equation is of the form: 

ap + bq = E{x,y). (m.2.8) 

Now let z = z(x, y) be a solution to: 

p = E{x,y) ("1-2-9) 
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-where x and y are defined by: 

x = a(x + y) (in.2.10) 

y = b[x — y). 

Then z(x,y) = z(x{x, y), y(x, y)) is a solution to (III.2.8). 
First we express x and y in terms of x and y: 

, x - ± + i Cm-2.il) 

2a 26 

V ~2a~ 26* 

Now we determine the first order partial derivatives of z with respect to x and y in 
terms of the first order partial derivatives of z; 

£! == 5?£f?£ + f?££^ = f?i2- + ^i- (IIL2.12) 

dx dx dx dy dx dx 2a dy 2a 

di_ tete,tedy == d£l__dz1_ 
dy ~~ dx dy dy dy ~~ dx 2b dy 26 

S + ^ = ?r = B(xry) (in.2.13) 

dx dy dx 
a |i + & |! = J B(x,y) (IH.2.14) 



Thus: 



dx dy 
which is the same equation as (HI.2.8). 

Case 2) a 7^ and 6 = 

The image irradiance equation is of the form: 

ap = E{x,y). (III.2.15) 

As a 7^ 0, we can write this equation equivalently as: 

p= E( x >y) (III.2.16) 

which is of the form (III.2.6). 



a 
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Case 3) a = and b j£ 

The image irradiance equation is of the form: 

bq = E{x,y). (III.2.17) 

As b ^ 0, we can write this equation equivalently as: 

q = E ( x > y y (m.2.18) 



b 

Now let z = z(x, y) be a solution to: 

P== —b— 
where x and y are defined by: 



(III.2.19) 



x = y (III.2.20) 

y = x. 

Then z{x, y) = z(x(y), y(x)) is a solution to (III.2.17). 

Now we show that the first order partial derivative of z with respect to y is the first 
order partial derivative of z with respect to x: 

dz = d 1 dx^ p (III221) 

dy dx dy 

Thus: 

dz_ = E{x,y) (ni 2 22) 

3y 6 

which is the same as equation (III.2.17). 

Hence it suffices to examine linear image irradiance equations of the form (III.2.6). | 

As previously mentioned, a FOPDE has, in general, infinitely many solutions. 
What kinds of conditions can one impose such that the solutions to a FOPDE are 
restricted to a single one? A solution to the Cauchy problem which is the problem of 
constructing an integral surface passing through any given curve C, provides us with 
one answer to this question and is stated in the following theorem: 

Theorem: Let p = E{x,y) a linear image irradiance equation and let C be an 
initial, continuous and non- characteristic curve. If A (III. 2. 3) does not vanish 
along C, then there exists a unique, smooth integral surface through C. 
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Proof: The proof follows from the existence and uniqueness theorem for ordinary 
differential equations and can be found in [COHI62b, pp. 145-147] and in section A.6. 
I 

If A vanishes along C, then C is either a characteristic curve and the equation has 
infinitely many solutions, or the integral surface does not have continuous derivatives 
along C. In particular at the end of section A.6 the following lemma is shown: 

Lemma: Let ap -f bq = E(x, y) be a linear FOPDE. Let C be a non- 
characteristic initial curve for which A vanishes (III.2.3). Then the solutions to 
this FOPDE are cylindrical surfaces perpendicular to the x-y plane. 

Proof: The proof can be found in section A.6. | 

In the previous theorem C was assumed to be a continuous curve. However, it can 
be shown that "any singularities of the initial data propagate in the x-y plane along 
the projection there of the relevant characteristic curve" [GAR64, p.22]. This is not 
surprising since characteristic curves can be viewed as branch curves along which two 
integral surfaces meet. 

Let x = x(i) and y = y(t) denote the parametric representation of a curve in 
the x-y plane. A point P = (x , yo) is called a double point, if several values of t 
correspond to P. In the case where an initial curve C or its projection onto the x-y 
plane has double points, the integral surface containing C has self-intersections and 
therefore z is not a single- valued function of i and y. 

The projection of a characteristic curve onto the x-y plane is called a base charac- 
teristic. In the case of a linear image irradiance equation the base characteristics are 
unique. An integral surface of a linear continuous image irradiance equation can have 
a bounding contour and we can show the following lemma: 

Lemma: Let ap -j- bq = E(x, y) be a continuous image irradiance equation and 
let an integral surface which has a bounding contour be defined by z = z(x, y). 
Then its b- silhouette is a base characteristic. 

Proof: For points On the bounding contour p and/or q are singular. Thus any initial 
curve which has a point in common with the bounding contour is discontinuous and 
the singularity propagates along a characteristic curve. | 

111.2.2. The General Continuous Image Irradiance Equation 

We discuss now general continuous image irradiance equations. A more detailed 
exposition of this material can be found in section A.3. In deriving the system of 
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characteristic equations (III.2.2) it is also assumed that an integral surface z has con- 
tinuous second order partial derivatives and therefore that the equation p y = q x holds. 
In [PLI54] it is shown that an integral surface can also be built from characteristic 
strips in the case when only its first order partial derivatives are continuous. 

The Cauchy problem for a general FOPDE is formulated in the same way as for 
linear FOPDE's and is discussed in more detail section A.7. Let C be an initial curve 
specified in parametric form by x = x{t),y = y{t) and z = z{t). Then p{t) and q(t) 
along C can be determined by solving the two equations: 

R{p(t),q(t)) = E(x(t),y{t)) (HI.2.23) 

As (III.2.23) is, in general, a nonlinear equation in p and q, several solutions may be 
possible for p{t) and q{t) along C. To "avoid inessential reference to possible multiple 
valuedness of solutions for p and q along C" [COHI62b, p.80] it is assumed that p and 
q are also known as initial data. For the following discussion we assume that p and 
q are specified along C. In other words the initial data is an initial strip (which we 
denote by d). Now if A 7^ (III.2.3) along C then a unique integral surface may be 
obtained: 

Theorem: Let R{p,q) = E{x,y) be a continuous image irradiance equation. 
Then for every initial strip which is not a characteristic strip and for which 
A^O, there exists a unique integral surface through this strip. 

Proof: The proof follows from the existence and uniqueness theorem for ordinary 
differential equations and can be found in [COHI62b, pp.79-82] and in section A.7. 
It has been shown [HA28] that if the initial data does not have continuous second order 
derivatives, then the solution does not have continuous first order derivatives. I 

An integral surface of an image irradiance equation is a solution to a Cauchy 
problem only under the assumption that R 2 p + R] ¥" (section A.3). Let S be an 
integral surface of an image irradiance equation for which the condition R\ + R q ^ 
holds. Then there exists an initial strip d such that the solution to the characteristic 
equations with C x as initial values is the surface S. Thus the points of a reflectance 
map at which R 2 p + R 2 q = have to be investigated separately. Such points are called 
stationary points and are defined as: 

Definition: A smooth function f{x,y) has a stationary point at (i ,J/o) if 

/*(xo,</o) = (HI.2.25) 

fy{XQ,yo) = 0. 
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We discuss now the conditions under which stationary points constrain the integral 
surfaces of an image irradiance equation. First we must define critical elements of the 
characteristic equations of an image irradiance equation: 

Definition: Let R{p,q) = E{x,y) be an image irradiance equation. Then a 

point denoted by {x Q ,y ,Po,qo) is a critical element if (x ,yo) is a stationary 

point of E(x, y) and (p 0> <7o) is a stationary point of R{p, q). 

The point (z , yd, Po, <?o) is a critical point of the image irradiance equation if it 

is a critical element and if the values (x , yo,Po, <?o) satisfy the image irradiance 

equation. 

The point {x ,y ,po,qo) is a singular point if it is a a critical point for which 

the values (po,9o) are uniquely determined by the values (x ,yo)- 

A point P is an isolated critical element if in some neighborhood of it, it is the 

only stationary point of E{x,y) and R{p,q). Isolated critical (singular) points 

can be defined similarly. 

Hence at a singular point the gradient of each integral surface can be uniquely deter- 
mined from the measured brightness. For example, the lines i = and p = consist 
entirely of critical points of the image irradiance equation p 2 = x . Note that these 
points are not singular as the value for q cannot be uniquely determined at a point 
(0,y) in the x-y plane. The point (x,y,p,q) = (0,0,0,0), however, is a singular point 
of the image irradiance equation p 2 + q 2 = * 2 + y 2 . Thus the tangent plane to each 
integral surface of this equation at the point (0,0, z) is parallel to the x-y plane. In 
chapter IV we will show how singular points constrain the integral surfaces of an image 

irradiance equation. rtDn „ . 

Discussing the Cauchy problem further, the case A = for a general FOPDfc, is 
analogous to the case A = for a linear FOPDE. If an initial strip is a characteristic 
strip, then the equation possesses infinitely many solutions. Again, the question arises 
as to whether it is possible to specify an initial strip C u (i.e., a curve C and p and 
q along it), such that A = and d is not a characteristic strip. Such a strip can 
be specified but "then there exists no integral surface which contains this initial strip 
and has continuous derivatives up to the second order in its neighborhood" [COHI62b, 
p 83] However, it might be possible to construct a surface through C x . Then the curve 
C is a singular curve for that surface [COHI62b, p.83], i.e., along this curve the function 
defining the surface does not have continuous second order partial derivatives. 

Thus the solutions to a continuous image irradiance equation do not necessarily 
have continuous partial derivatives everywhere. In particular p and q can be singular as 
shown in the example below. So, one integral surface of a continuous image irradiance 
equation can have a bounding contour although this does not imply that every integral 
surface of this equation has a bounding contour. For instance: 

z(x,y) = ix^ + f(y) (UI-2.26) 
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where / is any continuous function, defines an integral surface of the following image 
irradiance equation: 

(p 2 ~ I) 2 = (l-* ! ) 2 (III.2.27) 

(p' + l) 2 (I**!)'' 
Note that for x = 0, p, the first order partial derivative of z with respect to x as 
defined by (III.2.26), is singular. Thus this surface has a bounding contour and its 
image contains a b-silhouette. However E(x,y) and §f are continuous along the curve 

X = On the other hand, there exist integral surfaces of (III.2.27) which do not have a 
bounding contour. For example, the integral surface defined by: 

*(*,v) = f** + *(iO ( IIL2 - 28) 

where g is any continuous function is a solution to (III.2.27). 

There are, nevertheless, continuous image irradiance equations for which all integral 
surfaces have a bounding contour as shown in the next lemma. Any such equation can 
always be transformed into a singular image irradiance equation such that the original 
and the transformed equations have the same solutions: 

Lemma: Let R(p,q) = E{x,y) be a continuous image irradiance equation. 
Then all of its integral surfaces have a bounding contour, if there exist a G 
bijection g, a C 1 function / and if the following conditions hold: 

• The reflectance map can be written as: 

R(p, q ) = 9(f(p,<l))- ( IIU - 29 ) 

• There exist finite values fori and y denoted by x and y such that: 

lim 9 - 1 (E(x,y)) = ±oo. (HI.2.30) 

I-+XO 



x 
V-*VO 



Proof- The solutions to the continuous image irradiance equation R{p, q) — E{x, y) 
and the equation f(p,q) = f-Wz.y)) are the same. Note that the points (x ,yo) 
constitute the b-silhouette. Thus by virtue of the preceding discussion of bounding 
contours and b- silhouettes, the lemma is self-evident. I 

In section III.3 and chapter V we will discuss the solution to the reconstruction 
problem in the case where a b-silhouette can be uniquely identified in an image. 

The last case of interest here is where an initial curve C degenerates to a point P 
(see also section A.7 and note that the following discussion is relevant only for FOPDE s 
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which are not quasi- linear). Then "all characteristic curves through a fixed point P 
of x,2/,z-space form an integral surface" [COHI62b, p.83]. Such a surface S has a 
conical singularity at P and the Monge cone is the tangent cone to 5 at P. A surface 
constructed in such a manner is called the integral conoid of the partial differential 
equation. For a given image irradiance equation and a given point P there exists 
exactly one integral conoid. 

In summary, the method of characteristic curves can be used to determine the 
solutions to a FOPDE. An initial curve which lies on the integral surface must, in 
general, be specified in order to restrict the solutions to a continuous image irradiance 
equation to a single one. The critical points of an image irradiance equation have to 
be analyzed separately. 

111.2.3. Edges and Vertices 

In this section we give examples of how to formulate the Cauchy problem such 
that its solution is an integral surface containing an edge. As previously stated, if an 
initial curve C is discontinuous then the surface normals to the integral surface which 
contains C are discontinuous. We give a very simple example illustrating this. Two 
planes intersect along a straight line as depicted in figure 2. The image irradiance 
equation is assumed to be linear and the same for both planes: 

p + g = l. (IH-2.31) 

We pose the following questions: Do the planes have to be oriented in a particular 
way to assure that the previous image irradiance equation holds? Are there any other 
integral surfaces which contain the edge constituted by the intersection of two such 
planes? The answers to these questions depend upon additional information about the 
scene which is available. Referring to figure 2, if we specify the curve K in three space, 
it can be used as initial data to uniquely reconstruct the two planes from the image 
irradiance equation. As the image irradiance equation is linear, the curve T has to be 
a characteristic. We observed earlier (section III.2.1) that characteristic curves can be 
viewed as branch curves at which two different integral surfaces meet. Through every 
curve which is not a characteristic there exists exactly one smooth surface. Thus if 
an integral surface of (III.2.31) contains an invisible edge, its projection onto the x-y 
plane is a line whose slope is 1. (The characteristic curves of a linear image irradiance 
equation whose reflectance map is of the form R{p,q) = ap + bq where a and b are 
nonzero constants, are straight lines whose slope is proportional to £.) Recall that the 
linear reflectance map describes, for instance, the reflectivity properties of the maria 
of the moon and that the constants a and b specify then the direction towards the sun 
(section III.2.1). We conclude therefore that any ridges on the moon which are not 
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Figure 2. Two intersecting planes 



visible are in the direction of the sun. 

On the other hand, if only the edge T is specified, the two surfaces (which may 
not be planes) cannot be reconstructed in a unique way from a linear image irradiance 
equation. This follows from the fact that T has to be a characteristic curve. For 
example, the surfaces defined by the next two equations are solutions to (III.2.31): 

*(x,y) = x + (y-x) 2 (IH.2.32) 

z{x, y) = x + sin(y - x). (III.2.33) 



(III.2.34) 



At their intersection these two surfaces form an edge which is defined by: 

x = y z = x. 

Furthermore the two planes (which are also solutions to (III.2.31)) defined by: 

x _ 2 y + z = (IH.2.35) 

2x _ 3y + * = (III.2.36) 

give rise to the same edge. 

If an image irradiance equation is nonlinear, however, there is only a small number 
of surfaces which can give rise to a specific edge. In particular an edge cannot be a 
characteristic curve, as two integral surfaces which meet along a characteristic curve 
also have the same tangent plane there. So let T be an edge given in parametric form. 
Using equations (III.2.23) and (III.2.24) we can determine the possible values for p and 
q along T. Since the image irradiance equation is assumed to be nonlinear, several 
solutions for p and q are possible, so many different integral surfaces can be constructed 
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LINE1 



SURFACE 3 




Figure S. Vertex 



each of which has the edge embedded in it. Any two of these surfaces will intersect 
at an edge. Hence in the case of a nonlinear image irradiance equation, only a small 
number of different surfaces can give rise to a particular edge. 

Let us now specify a vertex V in space from which three straight lines emanate 
(referred to as a corner) as shown in figure 3. The equations of the three lines are 
given in three space. Assume, as well, that the following image irradiance equation 
holds in all three regions: 

p 2 + <y a = l. (III.2.37) 

We ask how many integral surfaces exist which have the same shading and contain the 
vertex and the three lines. (Note that these lines are not visible.) Clearly, all such 
surfaces will have a singularity at the vertex, i.e., the function z = z{x,y) defining a 
surface is not differentiable at V. A priori there are two ways of interpreting the lines 
leading from the vertex: 

1) as curves lying on a smooth surface or 

2) as edges. 



We will discuss both interpretations. 

Assuming case 1, we are faced with the problem of constructing a smooth (except 
at V) surface S which has a conical singularity at V and has the three lines embedded 
in it. Thus the vertex is a degenerate initial curve (see sections III.2.2 and A.7) and S 
is the integral conoid. If such a surface S exists, the image irradiance equation cannot 
be linear, the generators of the integral conoid must be straight lines and the vertex 
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cannot be a singular point of the image irradiance equation. In section III 2.2 it was 
shown that there cannot be another surface besides the integral conoid which has a 
conical singularity at V as the integral conoid is uniquely determined by the image 
irradiance equation and this vertex. So if the integral conoid with its vertex at V can 
be constructed such that the three lines lie on it, the surface S is unique. In the case 
where V is a singular point or the generators of the integral conoid are not straight 
lines, no surface S can be found and the three lines necessarily specify three edges. 

We show now that to be able to construct the surface S as described in the preceding 
paragraph, the three lines necessarily are characteristic curves of the image irradiance 
equation. So suppose that the integral conoid can be constructed at V such that it 
contains the three lines. Then the integral conoid and the Monge cone are identical, 
which follows from the fact that the Monge cone is the tangent cone to S at V. Ihus 
if the integral conoid has straight lines, they must coincide with the generators of the 
Monge cone. The lines leading from the vertex are therefore characteristic curves of 
the image irradiance equation. In other words, if the lines are characteristic curves and 
if V is not a singular point, then an integral conoid which contains the corner exists. 

In our particular example the integral conoid of (III.2.37) is a right circular cone 
whose generators are inclined to the x-y plane at 45 degrees. (A cone is called a right 
circular cone if the angle between its axis and a plane containing the circular cross 
section is 90 degrees.) Thus if the three lines are characteristic curves, the right circular 
cone is the desired surface S. 

The other interpretation of the corner is that the three lines originating at the 
vertex specify three edges. Again only in the case where the image irradiance equation 
is nonlinear can we possibly find three surfaces which form the three edges since, if 
the image irradiance equation is linear, different integral surfaces intersect only along 
characteristic curves. Yet the base characteristics are all parallel for a linear equation, 
hence the projection of the three lines onto the x-y plane cannot be base characteristics 
and so the lines cannot be characteristics. 

We want then to find three surfaces which intersect along the three lines using these 
lines as initial data. For a particular corner, three such surfaces do not nece^arily ^ exist 
as can be easily seen using equation (III.2.37). The integral surfaces of (III.2 37) that 
are constructed using a straight line (which is not a characteristic) as initial data are 
planes inclined at 45 degrees to the x-y plane. Although through three lines which 
constitute a corner we can always find three planes, they are not necessarily inclined 
at 45 degrees to the x-y plane and are therefore not necessarily integral surfaces of 
(III 2 37) However, if three integral surfaces do exist, they are unique. From the (two 
lines which are embedded in every surface, p and q can be uniquely determined. If 
these values for p and q as functions of i and y satisfy the image irradiance equation 
then each of the surfaces can be uniquely determined. 

Summarizing the observations made above: if three lines which form a corner are 
characteristic curves of an image irradiance equation, then the surface containing the 
corner is the integral conoid of the equation. Otherwise such a corner can be used to 
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reconstruct the three surfaces which form it. 

111.2.4. Gaussian Image 

In the previous sections we reviewed the well known results concerning continuous 
FOPDE's and showed that in general, a continuous FOPDE has infinitely many solu- 
tions. Only after imposing some additional constraints (in particular by specifying an 
initial strip which is not a characteristic strip) can the solutions be restricted to a single 
one. 

It is clear, then, that a continuous image irradiance equation does not contain 
sufficient information to uniquely reconstruct the shape of a surface from its image. 
What partial knowledge about a surface does an image irradiance equation provide? To 
understand this problem better, certain concepts taken from differential geometry are 
now introduced as they provide us with a formalism to discuss the shape of a surface. 
Some technical prerequisites are first reviewed briefly. A more detailed exposition can 
be found in any standard book on differential geometry, e.g., [CAR76]. For the rest of 
this section we assume that a function z = z{%, y) which defines a surface is C 2 . 

The Gaussian sphere is a sphere of unit radius. The surface normal at every point 
on the Gaussian sphere can have two possible orientations, referred to as the positive (or 
outward) and negative (or inward) directions. We adopt the convention that a surface 
normal on the Gaussian sphere points outward. The Gaussian mapping maps a point P 
on a surface into a corresponding point P G on the Gaussian sphere such that P and P G 
have the same surface normal. This mapping is well defined since no two points on the 
Gaussian sphere have the same normal. The Gaussian image (sometimes referred to as 
the spherical image) of a connected patch of a smooth surface maps into a connected 
region on the Gaussian sphere. The area of the Gaussian image of a surface is called 
its integral curvature. 

Points on a surface are classified according to the behavior of the tangent planes 
with respect to the surface: 

Definition: A point P on a surface is elliptic if the tangent plane at P does not 
intersect the surface at any other point. A point P is hyperbolic if the tangent 
plane at P intersects the surface in a curve which has two branches intersecting 
at P. A point P is parabolic if the tangent plane at P intersects the surface 
along a single curve (which can degenerate to a point). 



A surface which consists only of elliptic points is locally convex. We will say a 
surface is hyperbolic if it consists only of hyperbolic points. The points which separate 
regions where a surface consists of elliptic and hyperbolic points, respectively, are 
parabolic points. 

The grouping of points on a surface into elliptic, hyperbolic and parabolic points 
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Figure 4. Gaussian image of a neighborhood of an elliptic point 



can easily be expressed in terms of Gaussian curvature. Let F be a simply connected 
and bounded surface patch which is enclosed by the closed curve k and let G denote 
the Gaussian image of F, which is enclosed by the closed curve k G . "We divide the 
area G enclosed by k G on the sphere by the area F enclosed by k on the surface and 
then shrink the curve k down to a point P on the surface. Then F and G approach 
zero and their quotient approaches a definite limit K: 



lim ^ = K. 



(III.2.38) 



The number K defined in this way is called the Gaussian curvature" [HICV52, pp. 193- 
194]. A point P on a surface is elliptic if the Gaussian curvature is positive there, 
hyperbolic if the Gaussian curvature is negative, and parabolic if the Gaussian curvature 
is zero. If z = z(x,y) specifies a, surface and is C 2 , then the Gaussian curvature at a 
point (x, y, z) is defined by: 



K[x,y) = 



"zz-s-yj/ 



— Z 



xy 



(i + *2 + *!) a " 



(IH.2.39) 



The common notions of a surface being convex or concave do not refer to the type 
(as defined above) of points a surface consists of. They merely distinguish the two sides 
of a locally convex surface (or correspond to the two possible directions of a normal 
vector) with respect to a viewer. 

The mapping between a surface and the Gaussian sphere is one-to-one if the surface 
is either locally convex or hyperbolic. "If we move around an elliptic point along a 
small closed curve that lies on the surface, its spherical image - assuming that the 
surface has no double points - will also be a closed curve without double points, and 
this curve is traversed in the same sense as the original curve (figure 4). A small curve 
without double points about a hyperbolic point is also mapped into a curve without 
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Figure 5. Gaussian image of a neighborhood of a hyperbolic point 

double points, but in this case the sense is reversed (figure 5)» [HICV52, PP-^-l^. 
The spherical image of a surface which consists of elliptic hyperbolic and parabohc 
points consists of several sheets, i.e., several points on a surface get mapped into the 
same point on the Gaussian sphere. 

We will also need the definitions of a closed and a compact surface: 

Definition [CAR76, p.112]: Let A be a subset of 5ft 3 . We say that p <E » 3 is a 
limit point of A if every neighborhood of p in » 3 contains a point ******* 
from p. A is said to be closed if it contains all of its limit points A is bounded 
if it is contained in some ball of 5R 3 . If A is closed and bounded it is called a 
compact set. 

For example, the surface of a sphere is compact, whereas a paraboloid of revolution 
defined by z(x, y) = x 2 + y 2 is a closed but not compact surface. 

The notion similar captures formally what we mean by saying that two surfaces 
have the same shape: 

Definition: Two surfaces in JR 3 are similar if they can be mapped into each 
other by a composition of translations, rotations, reflections and dilations. 

111.2.5. Gradient Space 

In this section we will briefly discuss gradient space as popularized by Mackworth 
[MAC73i and Horn [H077] and which is now a standard tool in vision research. Our 
objective is to show how the concepts from differential geometry help us to understand 
the variety of integral surfaces of an image irradiance equation. We will also investigate 
how different constraints may restrict the number of possible solutions to a given image 
irradiance equation. 
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Figure 6. Map: reflectance map — ► Gaussian sphere 



Let 



R{p,q) = E{x t y) 



(III.2.40) 



be an image irradiance equation. For simplicity of exposition we assume that E(x, y) 
is denned at every point in the x-y plane. Then we can write the previous equation as 
a system of two equations: 



R(p,q) = c 

E{x,y) = c 



(III.2.41) 
(III.2.42) 



where c is a constant. In the p-q plane (also called gradient space), the graph of 
R(p, q) = c, for all possible values of c, is called the reflectance map. A reflectance 
map can be mapped onto the Gaussian sphere by placing the south pole of the Gaussian 
sphere onto the origin of the p-q plane by means of a simple projection from the center 
of the sphere, as illustrated in figure 6. Note that the mapping between gradient space 
and the Gaussian sphere is conformal. The graph of E[x, y) = c, for all possible 
values of c, can be drawn in the x-y plane, referred to as the image plane. For example 
consider the following eikonal equation: 



P a + * 2 = x a + y s 



(III.2.43) 



This equation can be rewritten as: 



c. 



(III.2.44) 



* 2 -fV 
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Figure 7. Graphs for the equation p a + q 2 — x* + y 



The graphs for these equations are depicted in figure 7. 

Let z = z(x, y) define an integral surface of an image irradiance equation. Then 
the gradient at every point of z is the vector (p, q, —1) where p = §* and q = §*. As 
discussed in chapter II, we assume that each point P = (x , y 0} z(x , y )) on the surface 
is mapped via orthographic projection into the point (x , y ) in the image plane. We also 
define a mapping from surface orientation to gradient space: The gradient (p , go, — 1) 
at P is mapped into the point (po, go) in gradient space. 

An image irradiance equation gives us some information about the correspondence 
between points in the image plane and points in gradient space. In particular, if (xo.yo) 
lies on the curve E(x,y) = c for some constant c , then (p ,9o) lies on the curve 
R(P, q) = co- Note that determining an integral surface (up to a constant factor) of 
an image irradiance equation is equivalent to specifying the correspondence between 
points in the image plane and points in gradient space. 

We want to show now that any two integral surfaces of an image irradiance equation 
need not be similar. For two such surfaces to be similar, it is necessary that they consist 
of the same type (i.e., elliptic, hyperbolic or parabolic) of points. An image irradiance 
equation, however, does not restrict the type of points on the integral surfaces. So 
by knowing one integral surface of an image irradiance equation, other such surfaces 
cannot, in general, be obtained through similarity transformations. To prove these 
assertions we examine some integral surfaces of (IH.2.43). In particular the surface 
defined by the following equation consists entirely of elliptic points: 



z(x,y) = 



* 2 + y 2 



(IH.2.45) 
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Figure 8. Graph for characteristic curves 



whereas the surface defined by: 



is hyperbolic. 



z = xy 



(III.2.46) 



As already noted several times, the basic tool for solving a FOPDE is the method 
of characteristic curves. Using gradient space, the characteristic curves of an image 
irradiance equation and the necessity of specifying an initial curve to restrict the possible 
solutions to a single one can be visualized easily [H075]. In the following we use 
only four of the five characteristic equations (III. 2. 2) of an image irradiance equation 
R{p, q) = E(x, y): 

dx 

= R p (p,q) (III.2.47) 



ds 

4 -«.(.,,) 



dq 
ds 



= E y {x,y). 



(III.2.48) 
(III.2.49) 
(III.2.50) 
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In the image plane, dx and dy can be viewed as the two components of a vector and 
E x and E y as the direction of a normal vector to the curve E(x,y) = c (where c is 
a constant). In the same way dp and dq can be interpreted in gradient space as the 
components of a vector and R p and R q as the direction of a normal vector to the iso- 
brightness curves R(p, q) = c. From equations (III.2.47) and (III.2.48) one can deduce 
that the two vectors (dx, dy) and (R p , R q ) are parallel. Similarly, equations (III.2.49) 
and (III.2.50) indicate that the two vectors {dp,dq) and (E x ,E y ) are parallel. 

Now let the point (x , J/o) m the image plane correspond to the point (p , <fo) in 
gradient space (see also figure 8 which is taken from [WOOD78, p.190]). Then a 
step in the image plane from the point (xo, t/o) iu the direction of a characteristic 
curve, corresponds to a step in gradient space from (po><7o) m * ne direction (E x ,E y ). 
Conversely, a step from the point (po>9o) hi the direction (dp, dq) corresponds to a 
movement in the image plane starting at (xq, j/o), in the direction (R p , R q ). Specifying 
an initial curve now gives a correspondence between a set of points, denoted by X, in 
the x-y plane and a set of points, denoted by P, in gradient space. If this initial curve 
is not a characteristic, the characteristic curves can be expanded from each point in X 
which just corresponds to expanding curves in gradient space from every point in P. In 
this manner a point (p, q) is assigned to each point in the image plane. Once p and q are 
known at every point on the surface, the function z =■ z(x, y) defining the surface can 
be found by integration (assuming that p and q satisfy the equation p y = q x ). Thus to 
determine a unique (up to translation in the z-direction) solution to an image irradiance 
equation, enough information has to be known so that the characteristic curves can be 
expanded simultaneously in the image plane and in gradient space. 

III. 3. The Singular Image Irradiance Equation 

In section III.l we classified image irradiance equations into two categories. In the 
previous sections, we discussed equations which fall into the first category, i.e., the 
case where E(x, y) is C 1 . This section deals with the second class, i.e., the case where 
E(x, y) is singular. 

Recall (chapter I) that a function E(x, y) is singular at a point (xq, j/o) if: 

Mm E(x, y) = ±oo. (IH.3.1) 

v-»»o 

Lemma: Let R(p,q) = E(x,y) be a singular image irradiance equation. Then 
the b-silhouette consists of those points (xoiVo) hi the x-y plane for which 
E(x, y) is singular. 



Proof: Each singular image irradiance equation defines a b-silhouette. The lemma 
follows. | 
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Let us recall here briefly the geometry of the image forming system as discussed 
in chapter II. It is assumed that the viewing direction coincides with the z-axis. 
Furthermore, we approximate the imaging situation by orthographic projection hence 
the lines connecting the viewer and points on the surface are parallel to each other and 
perpendicular to the x-y plane. For simplicity we will restrict our attention to images 
where a single b-silhouette occurs. Images with multiple b-silhouettes can be treated 
in a similar fashion. We will call a b-silhouette nondegenerate if it is a connected and 
smooth curve in the x-y plane. 

Let R(p,q) = E(x,y) be a singular image irradiance equation, let the b-silhouette 
be defined by w(x, y) = and let z = z(x, y) define an integral surface. Then 
its bounding contour consists of the points (x ,yo,z(xo,yo)) on the integral surface 
such that (x ,y Q ) belongs to the b-silhouette. We assume that integral surfaces are 
(piecewise) smooth (section III.l). So, no parts of an integral surface obscure each 
other and therefore the bounding contour is a set of (piecewise) continuous curves. As 
previously stated, the lines connecting the viewer and points on the bounding contour 
graze the surface. In other words, at every point on the bounding contour the tangent 
plane is perpendicular to the x-y plane. In terms of p and q this means that p and/or q 
assume infinite value at points on the bounding contour. The spherical image of points 
on the bounding contour (section III.2.4) corresponds to points on the equator of the 
Gaussian sphere. 

If the equation of the b-silhouette is known, the surface normal to points on the 
bounding contour can be determined. (This follows directly from the definition of 
bounding contour. The surface normal at a point (x ,yo,z(xo,yo)) on the bounding 
contour is {±Wx{zo,yo), zL w y{ x o>yo)}ty where the -f- and — sign distinguish the two 
sides of a surface with respect to the viewer.) This surface is tangent to a cylinder 
whose intersection with the x-y plane is the b-silhouette and whose generators are 
perpendicular to the x-y plane. 

So again we pose the question: What are the constraints necessary to obtain a 
unique solution to a singular image irradiance equation? We shall proceed in two ways. 
In section III. 3. 2 we show that the method of characteristic curves can be used to find 
the integral surfaces of a singular image irradiance equation. Unfortunately, specifying 
a bounding contour does not restrict the possible integral surfaces to a single one, as 
shown by means of an example in section III. 3. 2. In chapter V we investigate constraints 
which enable us to solve the reconstruction problem uniquely. 

The existence of a smooth integral surface is not guaranteed for every singular 
image irradiance equation as we now demonstrate. It is important to notice that we 
are looking for a global solution to an image irradiance equation, i.e., an integral surface 
which is defined at every point (x, y) for which the equation is defined. In particular, 
an integral surface should be bounded for points (x, y) which lie on the b-silhouette. An 
example of an image irradiance equation for which no global bounded solution exists 
is given by: 
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p=I. (IH.3.2) 

x 

The general solution to this equation is: 

z{x, y) = In x + f{y) + c (IH.3.3) 

where / is any continuous function and c is a constant. For points on the y-axis (i.e., 
x = 0) z(x, y) as defined by the previous equation is not bounded. 

For any given singular image irradiance equation, the points on the b- silhouette 
can be found by inspection of E{x, y). The following example should clarify this: 

«2_|_ 2 = g2 + y a — . (HI.3.4) 

P +9 l-(x*+y*) K ) 

Along the circle x 2 + y 2 = 1, E(x,y) assumes infinite value. An integral surface of the 
previous equation is a sphere defined by: 



z(x, y) = ^-(x' + y 2 ) (IH.3.5) 



with p and q given by: 



p = ~ x — (ffl.3.6) 

n/i - (* 2 + V 2 ) 
_ —y 

9 s/i-(x 2 -hy^)' 

Along the circle x 2 + y 2 = 1, p and q are infinite but the surface normal there is well 
defined as (x,y,0). Note also that if z = *(x,y) is a solution to an image irradiance 
equation, then so is z = z(x,y) + c for any constant c. As we are solely interested in 
determining the shape of surfaces, we will assume in the following chapters that c is 
given. (In other words, all integral surfaces are determined only up to translation in 
the ^-direction. Thus when we say, for instance, that there is a unique integral surface 
of an image irradiance equation we mean that it is unique up to translation in the 
z-direction.) In section V.2 we show that the convex and the concave hemispheres are 
the only two integral surfaces of (III.3.4) which have continuous second order partial 
derivatives. 

111.3.1. Marr's Occluding Contours 

One of the first to investigate b- silhouettes was David Marr [MA77]. His nomencla- 
ture differs from ours and to avoid confusion we compare the two terminologies. In his 
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Figure 9. Generalized cone 



paper bounding contours are referred to as contour generators, and b-silhouettes, as 
contours. 

Marr poses the question of what can be inferred about the shape of a surface 
from its b-silhouette alone. In order to more easily approach this problem, a priori 
assumptions about the surface must be made. The restrictions imposed are: 

Rl.) The surface is defined by a C 2 function. 

R2.) Each point on the bounding contour projects to a different point on the 
b-silhouette. 

R3.) Nearby points on the b-silhouette arise from nearby points on the bound- 
ing contour. 

R4.) The bounding contour is planar. 

As discussed in the paper, the first three restrictions are not severe. Marr points out 
that for bounded surfaces, R3 follows from R2 but he nevertheless chooses to state 
the two restrictions separately as "they have sufficiently different meaning" [MA77]. 
The third restriction states that there are no gaps in the viewing direction. Implied 
by this is that a bounding contour cannot be created by two surfaces, one of which 
partially occludes the other, as this would violate R2. A consequence of the first three 
restrictions is that a bounding contour is a continuous curve. 

The main technique used in [MA77] to interpret b-silhouettes is to examine their 
inflexion points, i.e., those points which separate convex regions on a b-silhouette from 
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concave ones. If a bounding contour is assumed to be planar, inflexion points on the 
b-silhouette imply the existence of inflexion points on the bounding contour. This 
observation leads to restriction R4 which as Marr observes, is a very strong one. He is 
therefore only able to analyze a small number of scenes. 

The main theorem of [MA77] (theorem 1) states that generalized cones are the only 
type of surfaces which satisfy the four restrictions above plus an additional assumption 
about the viewing direction. "A generalized cone (as shown in figure 9) is the surface 
generated by moving a smooth cross section p along a straight axis A. The cross section 
may vary smoothly in size (as prescribed by the axial scaling function h{z)), but its 
shape remains constant" [MA77]. Not only is R4 a very strong assumption, but for 
theorem 1 to hold, the viewing direction is confined to a plane which lies parallel to a 
cross section of the cone whose shape we wish to determine. However, a priori, there 
is no way to determine if the viewing direction satisfies the constraint of the theorem! 

The proof of theorem 2 as stated in [MA77] is wrong. This theorem claims that 
a surface obeys the four previously stated restrictions for ail viewing directions if and 
only if it is a quadratic surface. In the proof given it is assumed that the surface can 
be defined by a polynomial; the theorem is therefore only shown to hold for this special 
case. 

111.3.2. The Method of Characteristic Curves 

Let us recall here that the method of characteristic curves can be used to solve 
a continuous image irradiance equation. Is this method still valid in the case where 
an equation is singular? In fact, the answer is positive. The method of characteristic 
curves is based on a local theory as stated in [COHI62b, p.62]: "It should be emphasized 
again that all statements and derivations are in the small, i.e., they concern merely 
neighborhoods of points, etc., without necessarily specifying the extension of these 
neighborhoods." It is precisely the local nature of this theory which allows us to apply 
it to singular image irradiance equations. 

Let R{p,q) = E{x,y) be a fixed, singular image irradiance equation. The charac- 
teristic curves can be constructed in a fashion analogous to the continuous case. A 
difficulty arises only when the Cauchy problem is stated for this equation. In the con- 
tinuous case, for any strip denoted by C x (as defined in section III.2.2) which is not 
a characteristic strip and for which A ^ (III.3.3), an integral surface can be found 
which has Ci embedded in it. Furthermore this integral surface is uniquely specified by 
the FOPDE and the strip C\. The theorem used to obtain this uniqueness result is the 
existence and uniqueness theorem for ordinary differential equations. Unfortunately, 
in some neighborhood of a b-silhouette this theorem does not hold anymore, so there 
is no guarantee that the equation will have a solution for any initial conditions. 

To understand why the existence and uniqueness theorem for ordinary differential 
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(III.3.7) 



(III.3.8) 



Figure 10. Solutions to III.3.8 
equations does not hold at a singularity, we examine the equation: 

dy _ s + y 

dx x 

which is singular at (x,y) — (0,0). The set of solutions are: 

y = x ln(cx) 

where c denotes any constant. The graphs of some of the curves defined by the previous 
equation are depicted in figure 10 [SMI68a, p. 148]. 

The initial condition that the solution must pass through the point P, which has 
the coordinates (0, 0), does not suffice to pin down the solution to (III.3.7) uniquely. 
Furthermore there is no solution which passes through the point F° with the coordinates 
(0,3). This implies that for the initial condition P° no solution exists. In summary, 
one cannot claim that for any initial condition, (III.3.7) can be solved in a unique way. 
Moreover, the type of singularity in an ordinary differential equation constrains what 
initial data is consistent. 



If an image irradiance equation is singular, i.e., if E(x, y) is singular, then some 
or all of its characteristic equations are singular. Thus when specifying a bounding 
contour, there is no guarantee that a solution exists which has the given bounding 
contour embedded in it. Furthermore we can show the following lemma: 

Lemma: A singular image irradiance equation can have several integral surfaces 
each of which has the same bounding contour. 
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Proof: We will prove the lemma by showing it for a particular image irradiance equation. 
The two surfaces defined by (III.3.10) are solutions to (III.3.9) and each has the same 
bounding contour: 

p 2 + g 2 = — + 1 (m.3.9) 

it 
*(*,y) = -4**+y (m.3.10) 

z(x,y)= 4x* + y. 

I 

In chapter V we shall prove that the only integral surfaces of (III.3.4) are the 
convex and concave hemispheres. So for any bounding contour C other than the one 
specified by x 2 +y 2 = 1, z = 0, no integral surface exists which has C embedded in it. 
Note in addition that the concave and the convex hemispheres have the same bounding 
contour. By assuming that the surface is convex in the direction of the viewer, the 
solutions to (III.3.4) are restricted to a single one although this is not true if we forego 
this assumption, as shown in the proof of the previous lemma. 
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Chapter IV 



Singular Points 



IV. 1. Basic Concepts 

The question which we analyze in this chapter is: How much information about 
the shape of a surface can one obtain from a singular point of an eikonal equation? At 
a singular point of an image irradiance equation the gradient of all its integral surfaces 
is defined by the image intensities there (section III.2.2). As discussed in chapter I, 
equations of the form: 

p 2 + q*=E(x,y) (IV.1.1) 

are called eikonal equations and describe for instance the flux of the secondary electrons 
in a scanning electron microscope since this varies approximately as f(p 2 + q ) where 
/ is a continuous function [LAWH60]. 

To obtain our results we will impose some technical conditions upon E(x, y) (which 
we will discuss later in this section) and shall refer to eikonal equations which satisfy 
these conditions as constrained eikonal equations. The two results which we prove are: 

• There exist exactly two locally convex integral surfaces of a constrained 
eikonal equation in some neighborhood of a singular point. 

• At a singular point, the Gaussian curvature of each integral surface of a 
constrained eikonal equation has the same absolute value. 

The first statement can be expressed in other words as: if z = z(x, y) defines one 
locally convex solution, then z — —z(x,y) defines the other. Hence this result can be 
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viewed as a uniqueness result modulo the concave/convex ambiguity. The second result 
exhibits a situation where a "property" of surfaces can be determined from their given 
(common) image. Our theorems apply also to more general equations: in appendix II 
we exhibit a class of image irradiance equations whose solutions can be obtained by 
solving an appropriate eikonal equation. 

To characterize a constrained eikonal equation we need the following definition: 

Definition: A function g{x, y) vanishes precisely to second order at the point 
(0, 0) if its (limited) Taylor series expansion at (0, 0) is of the form: 

g(x, y) = ax 2 + /?xy + 7 y 2 + °((|*| + Ivl) 1 ) PV.1.2) 

where a, /? and 7 are constants, at least one of which is nonzero. 

Then a constrained eikonal equation is defined as: 

Definition: An eikonal equation p 2 + q 2 = E(x, y) is constrained if E{x, y) is 
a C 3 function satisfying the following conditions in some neighborhood of the 
point (io» Vq)' 

1) (xo.l/o) is a stationary point of E{x,y) 

2) E{x ,y o ) = (IV.1.3) 

3) E(x,y)>0 for (x, y) 7^ (x , y ) 

4) E{x, y) vanishes precisely to second order at (xo, yo). 

Let us discuss these conditions a bit further. The reflectance map of an eikonal 
equation is R(p, q) = p 2 -f- q 2 and therefore its stationary point is given by p = and 
q = 0. We have imposed the condition that {xo,yo) be a stationary point of E(x } y). 
Thus the point P = (x, y,p, q) = (x , yo, 0, 0) is a critical point of a constrained eikonal 
equation. Whence it follows from condition 2 that P is a singular point of such an 
equation. Prom the third condition we can deduce now that P is an isolated singular 
point (section III.2.2). By using a suitable linear transformation we may assume, 
without loss of generality, that the point (0, 0) is the stationary point of E(x, y). Since 
E{x, y) is assumed to be positive near the origin: 

az 2 + /?xy + <yy 2 >0 for (x,y)^(0,0) (IV.1.4) 

defines a positive bilinear form. Thus the subsequent inequality [BRSE75, p. 182] holds: 

a7 _^ >0 . (IV.1.5) 

4 
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Moreover, the constants a and 7 in the (limited) Taylor series expansion of E(x, y) must 
be positive. 

As mentioned earlier, if z = z(x, y) defines a locally convex solution, then so does 
z = — z(x, y). Hence, to simplify subsequent discussions, we will say that z(x, y) is 
a locally convex solution to a constrained eikonal equation if it satisfies the following 
positivity conditions in some neighborhood of the origin: 

1) 2(0,0) = 

2) z{x,y)eC 2 (IV.1.6) 

3) z(x,y)>0. 

Throughout this chapter it is assumed that a solution z = z(x, y) satisfies the positivity 
conditions and that the origin is the isolated singular point of a constrained eikonal 
equation. 

Horn [H075] used singular points to compute initial conditions sufficient to solve 
an image irradiance equation. In particular he assumed that an integral surface is 
convex at a singular point. In general, however, such a surface does not exist. Our 
results show when Horn's method can be used and we prove that in those cases, no 
initial conditions are needed to compute the convex surface. 

IV. 2. Preliminaries 

The results we prove in this chapter require a number of technical prerequisites 
which are introduced in this section. One of the key concepts is that of a Taylor series, 
which is discussed in order to deal with the problem of approximating a function by 
polynomials. To be able to write the Taylor series expansion of a C k function z{x, y) in 
a concise form we introduce the. notion of a homogeneous polynomial, which is a sum 
of terms of the same degree: 

Definition: A polynomial P(x, y) is a homogeneous polynomial of degree k if it 
is of the form: 

p(x,y)=j2 h i xj y k ~ j t™- 2 - 1 ) 

where for each j, hj denotes a constant. 

Definition: Let z(x, y) be a C k function. Then its (limited) Taylor series 
expansion at the origin is: 

k 
z(x, y) = z + J2 *i + <*(M + M)*) ( IV - 2 - 2 ) 
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where z denotes a constant and for each j, Zj is a homogeneous polynomial of 
degree j. For any j such that < j < k, there exists exactly one polynomial 
Zj which matches z at the origin up to the j-th derivative. 

Let us denote by / the first k + 1 terms in the Taylor series expansion of z(x, y) at 
the origin. Then / approximates z up to fc-th order, i.e., the error between / and z 
vanishes faster than any polynomial of degree k. 

Now we investigate the characteristic equations of a constrained eikonal equation 
and their solutions in some neighborhood of an isolated singular point. The four relevant 
characteristic equations are: 



(IV.2.3) 



dx 






~dl ' 


= 2p 




dy 

dt ' 


= 2q 




dp 
~dl ' 


= E x {x 


,v) 


dq 
~dl 


= E y {x 


,y)- 



dx 




~dl 


= 2p 


dy 




dt 


= 2q 



Since E(x,y) vanishes precisely to second order at (0,0), we can rewrite the charac- 
teristic equations in some neighborhood of the origin as: 



(IV.2.4) 

< ^ = 2ax + (3y + o{\x\ + \y\) 

^=/?x + 2 7 y + o(|*| + |2/|). 

One can view x, y, p and q as coordinates in 5R 4 . So, let £ denote the four- tuple (x, y, p, q) 
and let F(x, y, p, q) = p 2 + q 2 — E(x, y). We will say that £ 6 F if x, y, p and q satisfy 
the equation F = 0. Using this notation we introduce: 



Definition: Let p 2 -\- q 2 = E{x, y) be a constrained eikonal equation and let £ 
denote the four-tuple (x, y, p, q). The characteristic equations can be written 
as: 

f t =M + G[i) (IV.2.5) 
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where A is the four by four matrix: 



2 0> 
2 
A ~ ' 2a /? 



(IV.2.6) 



,p 27 o o; 

and where G has the following properites: 

1) G(£)<EC 2 

2) G(0) = (IV.2.7) 

3) f(0) = 

Every solution £ = £{*) to (IV.2.5) is called an orbit. The equation: 

^ = Ae ( IV - 2 - 8 ) 

is called the linearized characteristic equation. 

To describe orbits in some neighborhood of an isolated singular point the following 
concept is useful: 

Definition: Let 

at 
be an ordinary differential equation as in (IV.2.5). An orbit £ = £(*) is 
quasiradial if: 

lim e(0 = 0. OV.2.10) 

t-*±oo 

Equivalently we can express this definition as: 

lim (x(0,y(*)) = lim (p(t),q(t)) = (IV.2.11) 

t-*+oo t->+oo 

lim (s(t), y(*)) = lim (P(*).<7(0) = °- (IV.2.12) 

t-+ — oo t-» — oo 

In other words, if the characteristic curves are quasiradial, they are quasiradial in the 
image plane and in gradient space. Some quasiradial orbits are depicted in figure 11. 

The stable manifold theorem [ABMA80, p.527] and [HAR64, p.242] defines condi- 
tions under which the solutions to the linearized characteristic equations (1V.2.8) have 
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Figure 11. Quasiradial orbits 



the same topological structure as the solutions to the characteristic equations (IV.2.5). 
One of the crucial constraints is that all eigenvalues of A have non- vanishing real part. 
A precise definition of a manifold can be found in [ST069, p.203]. Here it suffices to 
observe that a surface defined by a C k function z = z(x, y) is a two-dimensional C k 
manifold. Before discussing the stable manifold theorem we state a theorem concern- 
ing two ordinary differential equations in two unknowns in some neighborhood of an 
isolated critical point: 



Theorem: (Node Theorem [HAR64, p.213]) Let 



dx 

It 
dy 

dt 



= a tl x + a 12 y + fi(x,y) 

= 021 X + 022V + /2O*, V) 



(IV.2.13) 



be a system of two ordinary differential equations where a,j for i,j = 1, 2 are 
constants such that 011032 — 012021 7^ and /, for i = 1, 2 are real continuous 
functions for which the following conditions hold in some neighborhood of the 
origin: 

h = o(l*l + \y\) h = o(M + |y|) (iv.2.14) 

Thus (x, y) = (0, 0) is an isolated critical point. If both eigenvalues of the 
linearized equation of (IV.2.13) are real and have the same sign, then all orbits 
are quasiradial and the critical point is called a node. In particular if both 
eigenvalues are negative, then each orbit tends to the origin as t — ► +00 and if 
both eigenvalues are positive, then each orbits tends to the origin as t — ♦ — 00. 
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The previous theorem can be viewed as a special case of the following theorem: 

Theorem (Stable Manifold Theorem): Let 

§=A£ + G(0 (IV-2.15) 

at 
be an ordinary differential equation where A is a constant matrix which has 
h eigenvalues with negative real part and d 2 eigenvalues with positive real 
part, where d and d 2 are both positive, and where G(£) satisfies the following 
conditions: 

1) G(£)GC 2 

2) G(0) = (IV- 2 - 16 ) 

3) f (0 = 0. 

Then there exist an e > such that (IV.2.15) has solutions £ = £(*) ^ 
satisfying: 

||£(*)||e £t -+ as t -+ oo (IV.2.17) 

and has solutions £ = £(t) ^ satisfying: 

||£(t)|| € £t -*■ as t-f —oo (IV.2.18) 

Furthermore, for sufficiently small e > 0, the point £ = and the set of points 
t(t) which satisfy (IV.2.17) sweep out a unique C 2 manifold of dimension dl, 
the stable manifold. Similarly, for sufficiently small e > 0, the point £ = and 
the set of points £{t) which satisfy (IV.2.18) sweep out a unique C manifold 
of dimension d2, the unstable manifold. 
Thus the curves which sweep out the stable (unstable) manifold are quasiradial. 

IV. 3. Integral Surfaces near a Singular Point 

We will now prove the main results of this chapter. 

Lemma: Let p 2 + q 2 = E{x, y) be a constrained eikonal equation. If a locally 
convex solution exists in some neighborhood of the singular point, then it is 
swept out by quasiradial characteristic curves. 

Proof: Suppose a locally convex solution z = z{x,y) exists. As z = z{x, y) is assumed 
to be C 2 , we can write p and q in some neighborhood of the singular point as: 

p = anz + any + o(|x| + |y|) ( IV - 3 - 1 ) 

q = ai2 x + a 22 y + o{\x\ + |y|) 
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where a tJ for i, j = 1, 2 are constants. As the origin is a singular point, p and q have 
no constant terms. Note that the Gaussian curvature K of z — z{x, y) at the origin is 
defined by: 

K = O11O22 — <*i2- (IV.3.2) 

Recall (section III.2.4) that a surface is locally convex at a point if its Gaussian curvature 
is positive there. Substituting the expressions (IV.3.1) into the first two characteristic 
equations (IV.2.3) of a constrained eikonal equation gives: 

% = 2( 0ll a: + a 12 y) + o(\x\ + \y\) (IV.3.3) 

at 

-^ = 2(a 12 x + a 22 y) + o(\x\ + |y|). 

Using the node theorem we deduce that the characteristic curves in the x-y plane are 
quasiradial if and only if both eigenvalues of the linearized equations have the same 
sign. A simple calculation shows that this is the case only when K > 0. Similarly, we 
can show that the characteristic curves in gradient space are quasiradial if and only if 
K > which in turn is true only if an and a 2 2 have the same sign. Assuming that 
K > 0, the sign of the eigenvalues is the same as the sign of an. | 

In the next section we will actually compute the coefficients a^ for i,j = 1, 2 such 
that K > 0. 

Finally we are ready to state and prove the main theorem of this chapter: 

Theorem: Let p 2 + q 2 = E(x,y) be a constrained eikonal equation. Then there 
exists a unique locally convex solution in some neighborhood of the singular 
point. 



Proof: It follows from the previous lemma that if a locally convex solution to a con- 
strained eikonal equation exists, it is swept out by quasiradial characteristic curves. So 
we have to show that such a solution exists and is unique. This is achieved by showing 
that the unstable manifold is the locally convex solution. To this end we investigate the 
linearized characteristic equations (IV.2.8) of a constrained eikonal equation. An easy 
calculation shows that the matrix A has two positive, real eigenvalues and two negative, 
real eigenvalues. Thus we can apply the stable manifold theorem, which states that 
there exist exactly two C 2 manifolds which are swept out by quasiradial characteristic 
curves. A locally convex solution therefore exists. From the node theorem we get that 
the solution z = z(x, y) satisfying the positivity condition is the unstable manifold, 
whereas the stable manifold is the surface defined by z = —z(x, y). Hence the locally 
convex solution is unique. | 
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IV. 4. Formal Power Series Solution 

In the previous section we-proved that there exists a unique, locally convex solution 
to a constrained eikonal equation in some neighborhood of the singular point, but did 
not show how to actually compute such a solution. We will do so in this section for the 
case where E(x, y) is analytic. In particular, we shall construct a formal power series 
solution. In the case where the eikonal equation is not analytic, such a solution tells 
us only about the behavior of an integral surface of the equation. 

Definition: A formal power series is an expression of the form: 

/ = *o+f>. (IV.4.1) 

where zq denotes a constant and for each j, Zj are homogeneous polynomials 
of degree j. We write the first order partial derivatives of / as: 

3=1 

where for each j, pj and qj are homogeneous polynomials of degree j. We will 
say that p 2 -f- q 2 = E(x, y) is a super constrained eikonal equation if E(x, y) is 
C°°. Then / is a formal power series solution to a super constrained eikonal 
equation if: 

(ir) 2 + (tt) 2 = ax2 + P*v + ™ 2 + I> ( IV - 4 - 3 ) 

dx dy *^ x 

where for each j, tj is a homogeneous polynomial of degree j. 

Suppose for a moment that we have computed a formal power series solution to a 
constrained eikonal equation. A theorem due to Borel [HAR, p. 261] states that there 
exists a C°° function z(x, y) which has a given power series as its Taylor series expansion. 
Unfortunately / does not determine z(x, y) uniquely since two different functions can 
have the same power series expansion. For example the functions g(x,y) and g(x, y) 
defined by (IV.4.4) have the same power series expansion at the origin: 

00 00 

g(x,y)=^ li + E^ ( IV - 4 - 4 ) 

3=1 j==i 

<i{x, y) = g{x, y) + e - ^" -f c - * 2 . 
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Note that the first order partial derivatives of g(x,y) and g{x,y) are the same only at 
the origin. Thus they cannot both be solutions to a given eikonal equation. As a formal 
power series solution satisfies the eikonal equation at the origin but not necessarily at 
any other point, it is not necessarily a solution to a given constrained eikonal equation. 
Yet the formal power series solution does tell us, qualitatively, about the behavior of a 
solution to an eikonal equation. Suppose a solution z to the eikonal equation exists in 
some neighborhood of the origin. Then z and 5 are tangent to at least second order at 
the origin. 

Lemma: Let p 2 -\-q 2 — E{x, y) be a constrained eikonal equation where E(x, y) 
is analytic. Then its formal power series solution is the solution to the equation. 



Proof: A version of the stable manifold theorem proves that if E(x, y) is analytic, then 
the stable (unstable) manifold is analytic [COLE55, p.330]. The lemma follows. | 

Lemma: Let p 2 + q 2 = E(x, y) be a super constrained eikonal equation. 
Then there exists a unique, locally convex formal power series solution to this 
equation in some neighborhood of the singular point. 

Proof: (outline) Equating the appropriate terms in (IV.4.3) we obtain 

• an equation for the quadratic terms and 

• a recurrence relation for each of the higher order terms. 

First (section IV.4.1) we shall prove that there is a unique solution to the equation for 
the quadratic terms if we impose the constraints that the formal power series solution 
be positive and convex. The next step (section IV.4.2) is to determine the higher order 
terms which is done by inductively solving the recurrence relation. If the quadratic 
terms have been determined such that the formal power series solution is convex, each 
step of this induction can be carried out uniquely. The first of the following two 
equations determines the quadratic terms and the second is the recurrence relation 
from which the higher order terms can be calculated: 

p 2 + g 2= a: r 2 +/?:ry + 7 y 2 (IV.4.5) 

2pip k + 2q x q k = ffk+i for fc > 1 (IV.4.6) 

where for each k, g k denotes a homogeneous polynomial of degree k. Each g k is easily 
computed using the power series expansion of E(x, y) as we will show in section IV.4.2. 
I 
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IV.4.1. Quadratic Terms 

In this section Aye determine pi and q\ for given values of a, /3 and 7 in (IV.4.5). 
Note that these terms define the Gaussian curvature of a solution to a constrained 
eikonal equation. Since p\ and qi are linear homogeneous polynomials they can be 
written as: 

pi = a n x + any (IV.4.7) 

9i = 021^ + 022!/ 

where the coefficients a,ij for i,j = 1,2 are constants. As px and q x satisfy equation 
(IV.4.5), the coefficients a tJ for i,j = 1,2 are constrained by: 

(o u x + a 12 y) 2 + (o 2 i x + o 22 J/) 2 = <*x 2 + 0xy + 7J/ 2 . (IV.4.8) 

Furthermore p\ and q x are the linear terms of the first order partial derivatives of a 
smooth function. Thus the so-called compatibility condition has to hold: 

d -l = |i (IV.4.9) 

dy ox 

which constrains the coefficients a\2 and a,2\: 

a 12 =a 21 . (IV.4.10) 

Equating appropriate terms in (IV.4.8) and using the compatibility condition we derive 
three equations for au,ai2 and 022 = 

°ii + °12 = <* 
2a 12 (a 11 +a 22 ) = /? (IV.4.11) 

a 12 +°22 — 7« 

This system of equations only has a solution if both a and 7 are not negative, which 
was assumed. 

As we wish to show that there exists only one positive convex formal power series 
solution z(x,y) to an eikonal equation, we want the coefficients a tJ for i,j = 1,2 to 
satisfy the previous equations and K to be positive: 

K = a n a 22 -a 2 12 . (IV.4.12) 

Thus the Gaussian curvature of z at the origin is determined only by p\ and q\. We will 
say that a solution for the coefficients a^ for i,j = 1, 2 is convex when K is positive. 
Note that for K to be positive it is necessary that an and 022 have the same sign. 
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Figure 12. Vectors p x and gi 



We can view pi and qi as vectors in the x-y plane. Then the vector product of the 
vectors p\ and q\ is: 



Pi X q\ — an<»22 — Hi- 
Recall that the vector product of two vectors can be written also as: 

Pi X?i = |pi|ki|sinp 



(IV.4.13) 



(IV.4.14) 



where <p denotes the angle between the vectors p x and qi as shown in figure 12. Thus 
for pi and q t to define a convex surface, the angle between them has to be less than 
180 degrees. We show now that the equations (IV.4.11) constrain this angle. The inner 
product of px and q\ is:' 

Pi • 9i = a 12 (an + o 22 ). (IV.4.15) 

Equivalently, we can write the inner product of pi and q\ as: 

Pi-9i = |Pilki| cos P (IV.4.16) 

where <p again denotes the angle between the vectors p x and q x . The following equations 
also hold: 



NHV^ii+^l 



(IV.4.17) 



|9l| = IV°12+ a 22l- 

So, we express the cosine of the angle <p between pi and q\ as: 



cosy? = 



2|v^7l 



(IV.4.18) 
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and there are two different angles <pi and <p 2 — 2* — fi which satisfy this equation. 
The next equation defines sin <p in terms of a 12 , a and 7: 

sin p = (±\/a - ^ 2 )(=bN/ 7 ~ af 2 ) ~ °L (rv . 4 . 19) 

Inasmuch as we are interested solely in convex solutions (i.e., <p < 180°), this equation 
has to be solved only for the case where the product of the square roots is positive, or 
equivalently, when an and a 22 have the same sign. Combining the last two equations 
yields: 

s/a-a\ 2 \/l-ah-a\ 2 = L JUL (IV.4.20) 

Since the the coefficients a, /? and 7 define a positive form, the previous equation is 
well defined and can be rewritten as: 

4( N /4a 7 -/3 2 + a + 7W2 = f- (IV.4.21) 

It is possible to determine a i2 in terms of a, /3 and 7 from this equation as long as 
a y£ 7 and (3 7^ 0. The case where = 7 and /? = will be investigated later. 

Case 1) a 7^ 7 and /? 7^ 

Let us denote the two possible solutions for a t2 obtained from (IV.2.21) by 012 and 
— ai 2 . Without loss of generality we assume that /? > 0. The results for /? < are 
analogous. The two solutions for an, 012 and a 22 are: 



on = \A* — a?2 °i2 = °12 °22 = Y 7 — °i2 (IV.4.22) 

a 11 = —\Ja — a\ 2 a 12 = — a 12 a 22 = — y7~ a i2- (IV.4.23) 

The coefficients defined by (IV.4.22) determine the unique positive convex solution. 

Case 2) a = 7 and /? = 

In this case equations (IV.4.11) can be written as: 

a\ x + a\ 2 = a (IV.4.24) 

2a 12 (a n +a 22 ) = (IV.4.25) 

a 2 2 + a* 2 = a. (IV.4.26) 
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Subtracting (IV.4.26) from (IV.4.24) we get: 

a?i-4»=<> flV.4.27) 

and deduce that: v 

a n =±a n . (IV.4.28) 

When an = —022. K is negative. So we only have to investigate the case where 
an = a 2 2- It follows from (IV.4.25) that a i2 = 0, thus we can express the coefficients 
a,j for i, y = 1, 2 in terms of a as: 

an = ±y/a 

ai 2 =0 (IV.4.29) 

Both solutions for the coefficients a*, for i,j = 1,2 are convex, but only one of them 
is positive. Thus we have shown that in the case where a = 7 and /? = a unique 
positive convex solution exists. | 

Actually we can also show the following theorem: 

Theorem: Let p 2 +g 2 = E(x, y) be a constrained eikonal equation. Then at the 
singular point, the Gaussian curvature of each integral surface has the same 
absolute value and is determined by the (limited) Taylor series expansion of 
E(x, y) at that singular point. 

Proof: Recall that the curvature at the origin, denoted by K, is: 

/f = ai 1 a 2 2-a? 2 . (IV.4.30) 

Using equations (IV.4.11) we derive an expression for this curvature in terms of a,/?, 7 
and aw- 







2 



(an -f- a J2 ) = - 
4a 



12 



(a„ + a 22 ) 2 = a - a 2 2 + 7 - 4, + 2ana 22 (IV.4.31) 

(an+a 22 ) 2 = a + 7 + 2K:. 



The desired expression for K is: 



K = L [ JL r -{a + 1 )). (IV.4.32) 

2 4af 



l 12 



55 

Prom equations (IV.4.18) and (IV.4.19) we get the two solutions for a\ 2 in terms of a,p 
and 7: 

P 2 (IV.4.33) 



°12 = 



„2 

°12 — 



4[ N /4a 7 -/5 2 + (a + 7)] 
P 2 



4[V4a 7 -/5 2 + (« + 7)]- 
Substituting these two expressions for a\ 2 into equation (IV.4.32) gives: 



2' 



Thus the absolute value of the curvature at the singular point of all integral surfaces of 
a constrained eikonal equation can be directly determined from the image intensities 
defined by E{x, y). I 

IV.4.2. Higher Order Terms 

In this section we calculate the higher order terms in a formal power series solution 
to a constrained eikonal equation using the solutions for Pl and q lt and show that th 
can be done in a unique fashion if Pl and 91 determine a locally convex solution. First 
we derive the recurrence relation which has to be solved to determine the higher order 
terms. Suppose that Pl and 91 are known. For p 2 and 92 the following equation must 

h ° ld: 2p 1 p 2 + 2« 7l < 72 = Z 3 . 0*4.35) 

For p 3 and 93 the following equation holds: 

2pip 3 + 2 9 i93 + V\ + q\ = U- (IV.4.36) 

Now, by assumption p 2 and q 2 are homogeneous polynomials of degree two. We can 
therefore write the previous equation as: 

2pip 3 + 29x93 = 94 (IV.4.37) 

where g 4 can be determined from h,p 2 and q 2 . Thus to determine, for each k, Pk and 
q k we have to solve the following recurrence relation: 

PiPk-i + 9i»-i = 9K (IV.4.38) 

where g k is a known homogeneous polynomial of degree fc. 

We now introduce a new coordinate system denoted by X and Y: 

X = Pl Y = qr. (IV.4.39) 
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We -can do so, as p± and qx are linearly independent vectors. (An easy calculation shows 
that if pi and q\ are linearly dependent they define a surface which has zero curvature 
at the origin. In other words the coefficients a,/? and 7 would not define a positive 
form.) Thus for each k > 2, we can write (IV.4.38) as: 

XP{X, Y) + YQ{X, Y) = F{X, Y) (IV.4.40) 

where F(X,Y) is a homogeneous polynomial of degree k and P(X,Y) and Q(X,Y) are 
homogeneous polynomials of degree k— 1. This follows immediately from the definition 
of X and Y. This last equation always has a solution, i.e.,: 

P(X,Y)= Fx{ ^' Y) (IV.4.41) 

since, using these two functions in (IV.4.40), we obtain: 

XF X {X, Y) + YF Y {X, Y) = kF(X, Y) (IV.4.42) 

which is Euler's homogenity relation [BRSE75, p. 246] for a homogeneous polynomial 
of degree k. 

Now we show that any solution for P{X,Y) and Q{X,Y) can be written as: 

P(X, Y) = P{X, Y) + C{X, Y)Y (IV.4.43) 

Q(X,Y) = Q(X,Y)-C(X,Y)X 

where C(X,Y) is a homogeneous polynomial of degree k — 2. Suppose that P 1 (X,Y), 
Q\X,Y) and P 2 {X,Y),Q 2 {X,Y) satisfy (IV.4.40). Then P[X,Y) = P 1 - P 2 and 
Q(X, Y) = Q 1 — Q 2 have to satisfy the following equation: 

XP{X, Y) + Y Q{X, Y) = 0. (IV.4.44) 

We are looking for two functions P{X,Y) and Q{X,Y) which satisfy the last equation 
and are homogeneous polynomials of degree k — 1. Furthermore P(X,Y] must vanish 
at Y = whereas Q(X,r) must vanish at X = 0. Thus P{X,Y) and Q(X,y) must 
be of the form: 

P{X,Y)= C{X,Y)Y (IV.4.45) 

Q(X,Y) = -C{X,Y)X 

which proves our assertion (IV.4.43). 
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Hence we can write the functions P[X,Y) and Q(X,Y), which satisfy equation 
(IV.4.40), as: 

p{XfY) = WOi + C (X,Y)Y (IV-4.46) 

Q(X,Y)= FY{X k >Y) -C(X,Y)X. 

Yet we want to find some P(X,Y) and Q(X,Y) which not only satisfy (IV.4.40) but 
also the compatibility condition: 

dP = dQ (IV.4.47) 

dy dx 

We will show that if Pl and q x define a convex solution, only one homogeneous poly- 
nomial C{X,Y) of degree k - 2 exists such that P and Q satisfy this condition. Using 
(IV.4.46), the partial derivatives of P and Q with respect to x and y are: 

to = k {Fxx ^y~ + XY W + *v l ' j + l ax a y + ay a y ' 
^ = l {FYX ~ax- + Fyy ax" )_ "aT C( ' ^ ( ax ax + ay ax > 

Recall the definitions of X and y in terms of x and y: 

X = a n x + a 12 y PV.4.49) 

y = ai 2 x + a 22y- 

So the partial derivatives of X and Y with respect to x and y are: 

SE aBu ^ = a 12 OV.4.50) 

ax ay 

ay dY 

ax ay 

As F(X,y) is a homogeneous polynomial, F XY = F rx . Using equations (IV.4.48) and 
(IV.4.50) we can express the compatibility condition (IV.4.47) as: 

-\F xx «i2 - Fyy*12 + F XY {a22 - an)] = (IV.4.51) 

- C(X, y)(a u + a 22 ) - C x {Xa n + ^12) - C Y {Xa 12 + roaa)- 

We wish to determine when there is only one polynomial C(X,Y) of degree k — 2 
which satisfies the previous equation. In so doing, we will compute the coefficients in 
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the polynomial C{X,Y). Then let 



k 



F(X,Y)='£hk-iX i Y k - i 

t=0 

k 

F X (X,Y)=J2fi>*-i iXi ~ lYk ~ i 

F Y (X, Y)=J^ f iik ^(k -yPY"-*- 1 (IV.4.52) 

F X x(x,Y) = £/,., fc _,i(i- i)x'- 2 y fc - 

*=2 

F y y (X, Y) = X; hk-i{k - i)(k - i - l)X«Y fc -*- 2 
»=o 

k— 1 

F X y(x,y) = J) /,■,»_,•(* -ipr'-^*-'- 1 . 

1=1 
Similarly, letting 

k— 2 

c(x,y)= ^•^^x'y*-'- 2 

J=0 
k 2 

C x (X,y) = X^.*-;-^- 1 ^ - '" 2 (IV.4.53) 

i=i 

c y (x, y) = X) c ;.*-;-2(* - i - 2pr'Y*-'- s 

3=0 



gives: 



fc— 2 

xc x (x,y)= X^^^yx'Y*-'- 2 

fc 2 

yC x (X, Y)=Y< Ci.k-j-dXi-'Y"-'- 1 (TV AM) 

XC y (X, Y) = X; Ciifc _ i _ 2 (A: - j - 2 )Xi+ l Y k -'-* 

3=0 
fc— 3 
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Using the last three systems of equations in (IV.4.51) we can rewrite the compatibility 
condition as: 

1 fc 

»'=2 
fc— 2 



£ f i>k _i{k - i){k - i - \)X i Y k - i ~ 2 a l2 + 

fc— l 

£ f itk -ii{k - i)X«- 1 r fc -'- 1 ( a22 - a n )] = (IV.4.55) 

t==i 

fc— 2 

- £ c j , k - j - a X*Y'-'-\a 11 + a 22 ) - 



fc— 2 

i=i 

fc— 2 

X;c il fc_ i _ 2 yx^- 1 y fc -^'- 1 a 12 - 

fc— 3 

£ Ci , fc _ i _ 2 (fc - j - 2^+ l r fc -'-»a ia - 
i=o 

fc— 3 

£ c;,fc-;- 2 (* - y - 2)X'Y fc " J - 2 a 22 . 
i=o 

Equating the coefficients of the powers X' x Y k ~' 1 for // = 0, ..., A; — 2 we obtain the 
following equations: 

For n = 0: 

-[/ 2)fc _ 2 2a 12 - f 0ik k{k - l)a 12 + /!,*_!(* - l)(a 22 - a n )] = (IV.4.56) 

— co.fc— 2(011 + 022) — ci.fc— 3012 — c ,fc_ 2 (fc — 2)a 22 

For < n < A: — 2: 

£[/,i + 2,*-m-2(/* + 2)(ju -f l)a 12 — / Ml fc_ M (fc — /x)(fc — m — l)ai2 + 
/m+i,*-m-i(* — M — !)(*: — M)( fl 22 — 011)] = (IV.4.57) 

— c /i,fc— m— 2(011 + a 22 ) — c Mifc _^_ 2 /xan — 

— c M +i,fc_ M _ 3 (M + l)ai2 — c M _i |fc _ M _i(A; — // — l)a a2 — 
c^.fc— M — 2 (fc — /j — 2)a 22 
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For /i = A: — 2: 

l[fk,ok{k - l)oi 2 - /fc_ 2f 22a 12 + /*-i,i(fc - l)(a 22 - a n )] = (IV.4.58) 

— Cfc— 2,o(Oll + O22) — Cfc_2,o(fc — 2)an — Cfc_3,i O12 . 

Denote by F M , for fi = 0, ..., fc — 2, the left hand side of each appropriate equation, 
multiplied by —1. Then we can rewrite (IV.4.56), (IV.4.57) and (IV.4.58) as: 

Fo =C0,fc— 2[ail + <»22(fc — 1)] + Ci,fc_3«12 

F* =c M -i,k_M-i(* - H ~ l)ai2 + (IV.4.59) 

C M ,fc_ M _2[ail(M + 1) + «22(fc — li — l)) + 
C M +l,fc— M — 3(M + 1)°12 
^fc— 2 =Cfc— 3,1«12 + Cfc_2,o[oil(fc — 1) + 022]' 

Now, for each fixed fc > 2 we need to determine the coefficients a t k—j—2 from the 
previous system of fc — 1 linear equations. This system can be written in matrix form 
as: 

A fc _!C = F (IV.4.60) 

and has a unique solution if the matrix Ak—i is nonsingular (i.e., its determinant is 
not equal to zero). Note that Ak—\ is a tridiagonal band matrix, so its determinant 
can be computed using the recurrence relation (IV.4.61) which we now derive. 

Let B n be an n X n matrix, let J9 n — 1 denote the submatrix of B n obtained from 
B n by eliminating its top row and leftmost column and let 2? n _2 be similarly obtained 
from B n —i. Then B n can be written as: 



(a b 0\ 
c £„-i J. 



(W.4.61) 



The determinant of B n may be computed in terms of the determinants of its sub- 
matrices: 

det B n = a det B n _i — be det B n -2 . (IV.4.62) 

It suffices to prove by induction that Ak—i is nonsingular if an 022 — «i2 > 0. 
Without loss of generality (see section IV.4.1) we can assume that an > and 022 > 0. 
The /xth row of Ak—i is: 

(k — ii — 1)012 au(l* + I) + *2i{k — fi — 1) ai2(/*-M)- (IV.4.63) 

Set j — k — fi — 2 and denote by A_, the determinant of the submatrix of Ak—i 
consisting of its j -j- 1 bottom rows and j -\- 1 rightmost columns. We prove that Ak—\ 
is nonsingular by induction on j. 
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For the basis case, we have to show that Ao 7^ and Ai 7^ 0: 

A = an(A:-l)+«2a (IV.4.64) 

A! = [aii(* — 2) + 2a 22 ][a n (fc — 1) + 022] — a 2 2 {k — 2). 

By assumption Ao is positive. Rewrite Ai as: 

Aj = a 2 n {k -2){k - 1) + 2a xx a 22 {k - 1) + 2<4 + (a xx a 22 - a 2 X2 ){k - 2). (IV.4.65) 

Since we have assumed that 011022 — a X2 > 0> & follows that Ai is always positive. 

Assume now that Aj_ 2 and Aj_i, for j > 1, are positive. We show that Aj is 
then also positive which completes the proof. Using (IV.4.62) we get: 

A, = [a xx {k - j - 1) + a 22 (j + 1)] A,_i - a{ 2 (k -j- l)jAj- 2 . (IV.4.66) 

where Aj_i is defined by: 

A y _! = [a xx {k - j) + a 2 2j]A J _ 2 - a\ 2 {k - j){j - l)A,-_ 8 (IV.4.67) 

and where we define A_i to be 1. Using (IV.4.67) in (IV.4.66) gives: 

A, =on(fc -j- l){[o u (fc - j) + a2 2 ;']A i _ 2 - a\ 2 {k — j){j - ljAy-a} + 

a 22 (j + l)Ay_i - a 2 12 (k - j - 1); Aj-2 (IV.4.68) 

Aj ={a xx a 22 — a\ 2 ){k — j— l)jAj- 2 + a 22 {j + l)A;-i + 

a xx {k -3- !)(*: - ;')l«ii A i-2 - «i 2 0' - l)Ai-s]. (IV.4.69) 

To verify that Aj is positive, one need only show, therefore, that: 

a u Aj-2 > a\ 2 {j - l)Ay_ 8 . (IV.4.70) 

We proceed, once again, by induction. For j = 3 we wish to show: 

onAi > 2a 2 2 A . (IV.4.71) 

Using (IV.4.64) and (IV.4.65) in (IV.4.71) this can be equivalently written as: 

anlaj^k - 2){k - 1) + 2a xx a 22 (k - 1) + 2a 2 2 + (a u a 2 2 - a\ 2 ){k - 2)] > (IV.2.72) 

■ 2a 2 2 [a xx {k — l) + a 22 ]. 

To see that the previous equation holds it suffices to observe that: 

2a 2 n a 22 {k — 1) > 2a 2 2 a xx {k - 1) (IV.4.73) 

2aua 22 > 2a 12 a 22 . 
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This solves the basis case. Assume, by induction, that the next inequality holds: 

auA i _ a >a; a (,--l)A i _ a . ( IV - 474 ) 

We need to show that: ,„. , „, 

anAi-i > <*? 2 JA,_2- f^- 4 - 75 ) 

Using (IV.4.67) in this last inequality gives: 
o2 1 (fc_j)A J _ 2 + ana2 2 ;A i - 2 -a 11 a2 2 (/ C -y)(y-l)A i _3 > a 2 12 jAj- 2 . (IV.4.76) 

This holds if: „ ,_ Tn __. 

aJiCife - y)A,_ 2 > o n aJ a (fc - j)(j - lJAy-s (IV- 2 - 77 ) 

which is true by the induction hypothesis (IV.4.74). 

Hence there exists a unique locally convex formal power series solution to a con- 
strained eikonal equation which can be easily computed. I 
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Chapter V 



B-Silhouettes 



V.1. Overview of Basic Concepts 

As discussed, our goal is to find necessary constraints such that an image can be 
interpreted in a unique way when its image irradiance equation is known. It was shown 
in chapter III that, in general, there can be an infinite number of different surfaces 
which satisfy the same image irradiance equation; in other words, the image of each 
of these surfaces is the same for a fixed imaging configuration. Recall that in sections 
III.2.2 and III.3 we discussed how to detect a b-silhouette in an image. It remains 
now to investigate whether the existence of such b- silhouettes can be used to interpret 
an image. In this chapter we will identify three constraints upon an image irradiance 
equation, one upon the reflectance map, one upon the b-silhouette and one upon the 
function E(x,y). If these constraints hold for some image irradiance equation, only one 
surface defined by a C 2 function which satisfies the equation exists. 

What kind of information about an integral surface can one deduce from a singular 
image irradiance equation? Let R{p,q) = E{x,y) be a fixed, singular image irradiance 
equation whose nondegenerate b-silhouette is defined by w{x, y) = 0. As previously 
discussed (section III.3), the surface normal to the bounding contour of an integral 
surface is parallel to {±w x , rbw v ,0). 

V.2. Uniqueness Theorem 

In this section we obtain constraints which assure that if an image irradiance 
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equation has a C 2 integral surface, it is unique. So let 

R(p,q) = E{x,y) (V.2.1) 

be a singular image irradiance equation. Consider the following constraints upon this 
equation: 

Cl) R(p,q)=P i +q 2 - 

C2) The b-silhouette denned by w{x, y) = is a closed, smooth curve in the x-y 
plane. Furthermore, the points {x,y) at which the image irradiance equation is 
defined, lie in the region hounded by this b-silhouette. 

C3) The function E{x, y) has exactly one stationary point (x , yo) and satisfies 
the following conditions in some neighborhood of (x ,yo): E(xo,yo) = 0, 
E{x,y) > for (x,y) ^ (x ,y ) and E{x,y) vanishes precisely to second order 

at(x ,yo). 

Uniqueness Theorem: Let R(p,q) = E{x,y) be an image irradiance equation 
for which constraints Cl, C2 and C3 hold and suppose a C 2 integral surface 
defined by z = z(x,y) of this equation exists. Then the only other solution to 
the equation is z = —z{x, y). 

Proof: Let R(p, q) = E{x,y) be a fixed image irradiance equation for which constraints 
Cl, C2 and C3 hold. First note that the point P = {x,y,p,q) = (x ,yo,0,0) is an 
isolated singular point of the image irradiance equation (see also section IV.2). There 
are then two observations which allow us to prove the theorem. First, as the b-silhouette 
is a closed curve, an integral surface of the equation has to be compact. Second, from 
constraints Cl and C3 we can deduce that such a surface is convex at the singular 
point, which allows us to apply results of the previous chapter. 

Suppose z = z{x,y) defines a C 2 integral surface of an image irradiance equation 
as defined in the uniqueness theorem. Then from C2 we may infer that z defines a 
compact surface (section ffl.2.4). Note also that z{x,y) is defined for every point (x,y) 
which lies within or on the b-silhouette and therefore has a bounding contour. Thus 
there exists a point P at which z has an extremum. In particular the tangent plane at 
P is parallel to the x-y plane. 

From condition C3 we can deduce that there exists a point {x ,yo,z ) on z sucn 
that the plane tangent to z at this point is parallel to the x-y plane, i.e., p(x ,yo) = 
and q(x Q , y ) = 0. By assumption this is the only singular point of the image irradiance 
equation and therefore the only point on z where the tangent plane is parallel to the 
x-y plane. Furthermore, as the image irradiance equation is singular, the point (x ,yo) 
lies in the interior of the region of the x-y plane bounded by the b-silhouette. Hence 
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P = (xo,yo,«o). By the assumption on E{x,y), z is either convex at P or has a saddle 
point. Since P is the point where the surface has maximal (minimal) height, z must be 

convex there. " . -• a m a r>t 

In chapter IV we proved that if an image irradiance equation satisfies Ol and ^6, 
there exists only one positive and exactly one negative convex solution denoted by z v 
and z n = —z p respectively. Thus there are exactly two integral surfaces z = z{x,y) 
and z = —z(x, y). I 

By using transformation methods, we can enlarge the class of singular image irradiance 
equations for which the uniqueness theorem holds. Let 

/( V + 2Bpq + Cq 2 + 2Dp + 2Eq) = E{x, y) (V.2.2) 

be a singular image irradiance equation Where / is a bijection and A,BC,D and E 
are real constants such that 6 > and AS < where 6, A and S are defined in the 
following equations: 

6 = AC — B 2 
A B D 



A= B C E 
DEO 

S = A + C 



(V.2.3) 



The constraints upon the constants A,B,C,D and E in equation (V.2.2) assure that 
the curves R(p,q) = c, for any constant c, are closed. Let the b-silhouette of the 
equation be a closed and smooth curve. Then (V.2.2) can be transformed into an image 
irradiance equation of the form (V.2.1) for which C2, holds as is shown in appendix 
II. If, after the transformation, the function E{x,y) satisfies C3, then the uniqueness 
theorem holds for (V.2.2). 

The next corolla ry follow s directly from the uniqueness theorem in this chapter. 
We will abbreviate \Jx* + y 2 by r. 

Corollary: Let p 2 + q 2 = E(r) be an image irradiance equation where E{r) 
satisfies constraints C2 and C3. Suppose a C 2 integral surface of this image ir- 
radiance equation exists. Then it is wtationally symmetric and can be obtained 
by integrating E{r). In this case the b-silhouette is a circle. 

Proof: First we write the eikonal equation in polar coordinates: 

4 + ?4 = m (V.2-4) 
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Let z = z(r) define the rotationally symmetric integral surface of the above eikonal 
equation. Thus z#(r) = and we can compute both rotationally symmetric solutions by 
integrating ±\/®W- It follows from the uniqueness theorem that the image irradiance 
equation has only rotationally symmetric integral surfaces. | 

Note that the above corollary does not hold if the image does not contain a 
b-silhouette. In section III.2.5 we showed that the integral surfaces of a continuous 
rotationally symmetric eikonal equation are not themselves necessarily rotationally 
symmetric. The uniqueness result for the special image irradiance equation (III.3.4) 
has been independently obtained by [DSY80]. 

V.3. Counterexamples 

In the previous section we discussed sufficient constraints under which the solution 
to a singular image irradiance equation is unique. Are these constraints necessary? 
Although we are not able to answer this question completely, we now shed some light 
upon it. In particular we try to find the class of image irradiance equations for which 
most likely there is no set of constraints that assure the existence of only one global 
solution. 

Image irradiance equations satisfying the constraints of the uniqueness theorem 
have closed iso-brightness curves, i.e., the curves R(p,q) = c are closed. So let us 
examine singular image irradiance equations whose iso-brightness curves are not closed. 
Such an image irradiance equation is given by: 

r + q = -(* + y)_ , (V.3.1) 

While constraint C2 holds for (V.3.1), an image irradiance equation where the reflectance 
map is R{p,q) = p -\- q never has a singular point. The general solution to (V.3.1) is: 



z(x, y) = ^/l-^ + y 2 ) + w{y - x) (V.3.2) 

where w is any C 1 function. Figures 13 through 16 illustrate some solutions to equation 
(IV.3.1). 

Another example of an image irradiance equation whose iso-brightness curves are 
not closed is given by: 

"-infer ™ 

This equation satisfies constraint C2, i.e., its b-silhouette is a closed and smooth curve. 
Furthermore (0,0) is the singular point of (V.3.3) and E(x,y) vanishes precisely to 
second order there although E(x, y) is not positive in the neighborhood of (0,0). One 
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Figure 13. z{x, y) = y/l — {x 2 + y 2 ) 




Figure 14. z(x, y) = y/l - (i 2 + V 2 ) + (y — *) 
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Figure 15. z{x, y) = y/l - (x 2 + y 2 ) + (y - x) 2 




Figure 16. z{x, y) = y/l - {x 2 + y 2 ) + (v - a:) 3 
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of the solutions to (V.3.3) is the sphere: 



^^l-^H-V) (V-3-4) 



whereas another surface satisfying (V.3.3) is: 



z{x,y) = f{t) + x 2 -y 2 where ( V - 3 - 5 ) 

t = 1 — (i 2 + J/ 2 ) a nd 



Recall that constraint C2 expressed the fact that the b-silhouette is a closed and 
smooth curve. We now demonstrate that if the b-silhouette does not obey C2, our 
uniqueness result does not hold. Previously, it was not expected that the , uniqueness 
results for image irradiance equations containing closed b-silhouettes would be different 
from the results for those containing open b-silhouettes. An example of an equation 
for which CI holds, but whose b-silhouette is not a closed curve is: 



*+o"«£ + l. ( V - 3 - 6 ) 



Equation (V.3.6) does not have a singular point. A solution to this equation is: 



z (i > y) = >A + y (v - 3,7) 
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Figure 17. z{x, y) = \/x-\-y 
which is shown in figure 17. Two other solutions to (V.3.6) are: 

*(*, y) = *(/£ + !) + I N^+i + 1) - m7~ + i - !)]• ( y - 3 - 8 ) 
*(*,y) = 



2 4 ^ VS* 



A more complicated counterexample is: 

p 2 + q' = 



2-u„2 = * 2 + y 2 , 

l — xy 



(V.3.9) 



This equation clearly satisfies CI but its b- silhouette is not a closed curve. The function 
E(x, y) vanishes precisely to second order at the singular point of (V.3.9) and is positive 
in the neighborhood of the origin. Two different solutions to (V.3.9) are given by 
equation (V.3.10) (shown in figure 18) and (V.3.11) (shown in figure 19). 



z{x,y) = 2y/l — xy 



z{x,y) = {l — xy) ] 



l — xy 



— 1 — atani 



1 x 2 — y 2 
— H 7T- 



1 — xy 



(V.3.10) 
(V.3.11) 
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Figure 18. z{x, y) = 2^1 — xy 



Note that the surface defined by (V.3.11) is C* 1 at the origin (although not C ) and is 
only defined in two quadrants, whereas the image irradiance equation is defined in all 
four quadrants. The surface defined by equation (V.3.10) is not convex at the origin. 
Only when the b-silhouette is a closed curve can we deduce that a surface is convex 
at the singular point, an observation which allows us to prove the uniqueness theorem. 
For the following eikonal equation we give two solutions which are both C 2 at the origin 
but only one of which is convex: 



p 2 + </ 2 = 



x 2 + y 2 
l — iV' 



(V.3.12) 



Equation (V.3.12) has a singular point and E{x, y) satisfies constraint C3. Two solutions 
to this equation are: 



z(x, y) = arcsin(iy) 



(V.3.13) 



z(x, y) = \/l — x 2 y 2 + 



* 2 -y 2 
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Figure 19. z{x, y) = (1 - xy)^/^ - 1 - atan^/j^ -T + 



2 a 
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Chapter VI 



Numerical Methods 



VI. 1. General Remarks 

In this chapter we discuss some numerical methods proposed to solve the shape 
from shading problem. The first shape from shading algorithm was implemented by 
Horn [HO70], who proposed to solve the characteristic equations (HI. 2. 2) using standard 
numerical methods for solving ordinary differential equations. Horn observes that the 
surface gradient at a singular point P (section III. 2. 2 and chapter IV) is uniquely defined 
by the image irradiance equation. By assuming that a surface is locally convex in some 
neighborhood of P and then estimating the curvature of this surface, he is able to 
calculate an initial curve (sections A.6 and A.7). Horn notes that using this heuristic, 
only one surface is calculated. However, it is not guaranteed that there exists a solution 
to any given image irradiance equation which is locally convex in some neighborhood of 
P and so the surface determined using Horn's algorithm may not be an integral surface. 
Furthermore, the algorithm computes at most one of the possible integral surfaces of 
an image irradiance equation. We proved in chapter IV that in the case of an eikonal 
equation where some technical conditions are imposed on E(x, y) a unique locally convex 
solution exists in some neighborhood of a singular point. In this situation then, Horn's 
algorithm computes the unique convex solution. (Note that using our results it is not 
necessary to estimate the curvature at a singular point in order to uniquely compute 
the surface.) Unfortunately, the algorithm is slow, numerically unstable, and relies on 
the presence of singular points. 
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• Strat [STR79] developed an iterative shape from shading algorithm which we 
discuss in detail in section VI.2. One of the major shortcomings of this algorithm 
is that an initial curve lying on a surface is required as an input. 

Brooks [BR079] suggests a Waltz-like [WA75] constraint propagation scheme to 
determine the shape of a surface from its shading. Like Strat, Brooks assumes that a 
surface is C 2 and that an initial curve which is embedded in the surface is known. 
In addition, he imposes an upper limit upon the curvature of a surface that can 
be determined. This constraint (which stems from experiments with human visual 
systems) is needed for the algorithm to converge. 

The three methods mentioned above can only solve the shape from shading problem 
under the assumption that a surface is smooth and does not have a bounding contour. 
If an image contains a b- silhouette it can be analyzed by using the algorithm due to 
Ikeuchi and Horn [IKH081] which is explained in section VI.3. 



VI.2. Strat's Algorithm 

Strat [STR79] proposes an iterative algorithm to solve an image irradiance equation 
under the assumption that an integral surface is defined by a C 2 function z = z(x, y). 
He assumes that image irradiance is measured at discrete points on the image plane 
and that the reflectance map is given. In this scheme, a square grid is imposed on the 
x-y plane. We denote a grid point by the tuple (i, j), the image irradiance measured at 
(i,j) by Eij, the first order partial derivatives of z = z(x, y) (which defines an integral 
surface) at a point (i,j) by p,- (J - and qij, and the reflectance map evaluated at (Pi,j,Q%,j) 
by Rij. The objective of the algorithm is to calculate p,- (J - and qij at every point (i,j). 
Strat does not address the problem of how to compute the function which defines the 
surface, i.e., Zij, from the values of p t) y and q^j. 

So let R(p, q) — E(x, y) be a given image irradiance equation. We now show how 
to obtain a function e (called the error function) in the variables p t|J and g tiJ which is 
zero when the image irradiance equation is satisfied at every point and pij and qij are 
the first order partial derivatives of z = z(x, y). Strat's algorithm computes iteratively 
the values for p tiJ and q^j which minimize the function e. This function is the weighted 
sum (over all (i, j)) of two functions e£ ■ and e* • whose derivations we now explain. The 
function ej •, defined by (VI.2.1), is zero when pij and g t(J satisfy the image irradiance 
equation: 

*ij = &J-Rij) 2 - (VI.2.1) 

The values for p, (3 and q it j which minimize e[ ■ cannot, however, be chosen inde- 
pendently; they are the first order partial derivatives of a C 2 function z = z(x, y). 
Consequently, equation (VI. 2. 2) holds for p and q [DIR72]: 



p y = q x . (VI.2.2) 
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Figure 20. Template 



Stoat assumes that the image irradiance equation is solved over a simply connected 
region in the x-y plane. Hence (VI.2.2) can also be expressed in integral form piR72|: 



/ 



pdx -f- qdy = 0. 



(VI.2.3) 



We now show how a discrete approximation of equation (VI.2.3) is used to derive e^. 
Strat suggests the loop around four adjacent points on the grid as a path over which 
the integral in (VI.2.3) should be taken. Such a path is depicted in figure 20. A discrete, 
first order approximation to equation (VI.2.3) is then given by: 

p itj + p i+ltj + qt+ij + qi+i,j+i - Pi+ij+i - PU+i - *.i+ 1 ~ qi >> = °- (VI-2 - 4) 
Thus the values for p x 3 - and q itj which the algorithm calculates must not only satisfy 
the image irradiance equation, but also equation (VI.2.4). In general, a point (t,jj 
belongs to four templates. For each template, a discrete approximation of (VI.2.3) can 
be derived and p ifi and q iti must satisfy each such approximation. So e^ is defined 
as: 

=(p,-+i,j + 9i+i,i + tt+i,j+i ~~ P»+U+i ~~ 

Pi,j+i — 9i,i+i + P«.i "~ 9*,i) + 

{pij-i + Pi-fij-i + q%+i,j—i + 9.+1.J — 

Pi+i.i — 9*,j— l — PiJ — 9t,i) + 

(pi_ li: ,_i + p it j-i + 9,,i-i — Pi-U — 

qi-\,j — 9«— i,i— i — P*-J + *.i) + 
(Pi-i,i + 9io+i — P»,i+i — Pt-i,i+i — 
9i— l.i+i — 9t-i,i + P«.i + 9*.i) • 



«.j 



(VI.2.5) 
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The total (global) error e is then defined as: 

c = XX+Ki) (VI.2.6) 



*,3 



where p is a constant scalar chosen appropriately to equate the dimensions and to 
weigh the magnitudes of the errors e^ and ej^. At each iteration, Strat's algorithm 
minimizes e. Hence at every point (i, j) the following two equations must be solved for 
Pi,j and q itj : 



de 

d Pi,i 
de 



= (VI.2.7) 

= 0. 



The Gauss Seidel method is used to solve equations (VI.2.7). 

It remains to describe how the initial assignment for the values ptj and qij is 
chosen. Strat assumes that the algorithm is applied to solve an image irradiance 
equation over a rectangular region in the x-y plane. For points (t, j) on the boundary 
of the rectangle, the values for p t>J and g t(J are given as input data and for interior 
points, pi t} and q it3 - are set to zero. As mentioned in section VI. 1, Strat's method is not 
very useful in practice, since it requires, as input, the initial values for p, (J and g t>J at 
the boundary of the domain in which the algorithm is applied. It is not shown whether 
the solution computed by the algorithm is independent of this initial choice of values 
for pij and qij at the interior points. 

Experimental results appear to support the conjecture that Strat's algorithm con- 
verges, although a proof of convergence is not given. One of the shortcomings of Strat's 
algorithm is that its performance depends upon the order in which the grid points are 
scanned, i.e., the order in which the values for p, (J and qij are updated. 

VI.3. Ikeuchj and Horn's Algorithm 

Ikeuchi and Horn [IKH081] designed and implemented an iterative algorithm which 
differs in various respects from the ones described in the previous sections. 

• Their assumption about the surface whose shading is to be analyzed is weaker. 
They observe that for a surface to look smooth it suffices for the function 
defining it to be only C 1 . Recall that under this assumption, an integral surface 
of an image irradiance equation can be built from characteristic strips (see 
chapter III and appendix I). 

• They show that a b-silhouette can be used as initial data, which is not possible 
using Brooks' or Strat's algorithm. 
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One of the strongest features of the method is that it incorporates data obtained from 
a b-silhouette. Recall that the surface normal at a point on the bounding contour is 
parallel to the normal vector to the b-silhouette. Before proceeding to describe and 
analyze Ikeuchi and Horn's algorithm, we review their slightly different formulation of 
an image irradiance equation. 

In chapter III it was shown that an image irradiance equation defines a relation 
between the image plane and gradient space. In the original formulation of this equation 
[H075] the (p,q) coordinate system was chosen to specify gradient space. With this 
underlying coordinate system an image irradiance equation can be formulated as a 
FOPDE. On the other hand this system has the disadvantage that for points on the 
bounding contour p and/or q assume infinite value. In order to find a numerical solution 
to an image irradiance equation it would be useful to transform gradient space into 
a different space SP, where (u,v) is the underlying coordinate system, so that the 
components of a surface normal vector (which always have finite value), instead of p 
and q, can be computed. To achieve this we have to find a space SP and a mapping 
m between gradient space and SP, satisfying the following constraints: 

• The mapping m between gradient space and a region of SP is one-to-one and 
onto. 

• The mapping m maps the whole (unbounded) p-q plane into a bounded region 
in SP. 

Recall (section III.2.5) that gradient space can be obtained by projection of the Gaussian 
hemisphere from its center onto a plane placed at its south pole. Ikeuchi and Horn 
observe that by choosing a different projection of the Gaussian hemisphere onto a plane, 
each point on this hemisphere gets mapped into a point other than infinity on the plane. 
Several projections are suggested, one of them being the stereographic projection. We 
can think of this projection in geometric terms as a projection onto a plane tangent to 
the hemisphere at the south pole. The center of projection is the north pole, not the 
center of the sphere. In this case the resulting coordinate system, denoted by (f,g), is 
defined in terms of p and q by: 

■ frfyi+p' + «*-!) (VL3.X) 

2<7(v/l+p 2 -fg 2 -l) (VL32) 

We can now write an image irradiance equation in the form: 

R(f,g) = E(x,y). (VI.3.3) 

Next we describe Ikeuchi and Horn's method of solving the shape from shading 
problem. The basic concepts of this algorithm are similar to those described in sec- 
tion VI.2. The inputs to this algorithm are the measured image irradiance and the 
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reflectance map R{f,g) and the desired outputs are the values f itj and g itj . (The 
problem of how to determine z ifj from the values of f itj and g ifj is not addressed by 
Ikeuchi and Horn.) There are two major differences between this and Strat's procedure. 
First, a different error function is used, and second the initial values for f itj and g itj 
are computed for points on the b-silhouette. (The case where there are no b-silhouettes 
in the image is discussed later.) For points not on the b-silhouette, f if j and gi,j are 
initialized to zero just as in Strat's algorithm. 

We now explain Ikeuchi and Horn's derivation of the error function. They exploit 
the assumption that the function z = z{x,y) defining an integral surface is C . A 
function F = F{x,y) is continuous at a point (x ,y ) if, for any e > 0, there exists a 6 
such that, for any (x, y) which lies in the circle defined by (x — x ) 2 + (y — y ) < *> 

we have* 

|F(x,y)-F(z ,yo)l<e. C^ 1 - 3 - 4 ) 

A discrete approximation of equation (VI.3.4) is given by: 
The error function e is therefore defined to be: 

c =D<i + Ki) where (VL3,6) 

jjHXj-Bur 2 (VI - 3J) 

(fc+i j + w-i j + ft j+i + ft.i-i ~ ^w)*- (VI,3,8) 

In equation (VI.3.6) X denotes a scalar factor chosen appropriately to equate the 
dimensions and to weigh the magnitudes of the errors e\ ti and t\ tj . 

As in Strat's algorithm, the function e is minimized. Thus, at every iteration the 
next two equations are solved for /,,/ and gi t f. 

de =0 (VL3.9) 



de 



= 0. 



d 9i,j 

Once again the Gauss Seidel method is used to solve equations (VI.3.9). A lower bound 
on the number of iterations needed using this method is proportional to the mesh size 
but an upper bound is not known. However, numerical experiments indicate that indeed 
the algorithm converges. 

As mentioned before, f itj and g itj are initialized to zero at every point except those 
on the b-silhouette. It is not shown that the algorithm is independent of these initial 

values. 

In the case where there are no b-silhouettes in the image, heuristics have to be 
used to initialize the algorithm. Ikeuchi and Horn do not specify the number of points 
at which the exact values for /,-,,- and g itj need to be known in order to guarantee that 
the algorithm will compute an answer. 
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Chapter VII 



Conclusion 



In this report, we have investigated the question of how much information con- 
cerning the shape of an object can be deduced from its shaded image. Even assuming 
that adequate data is available to derive an image irradiance equation is insufficient to 
solve the reconstruction problem uniquely; in general, for a fixed imaging configuration 
there are many surfaces which have the same shaded image. Thus our goal has been 
to identify constraints by which the reconstruction problem can be solved uniquely. 

We first analyzed the continuous image irradiance equation. The information 
needed to restrict the solutions to such an equation to a single one, was previously 
known [COHI62b]: If a strip is specified which is not a characteristic strip and along 
which A 7^ (III. 2. 3), then there exists a unique surface which contains this strip 
and whose image has a particular shading for a fixed imaging configuration. We were 
interested in whether edges could constitute an initial curve. In the case where the 
image irradiance equation is linear, it is not possible to distinguish among the multiple 
surfaces which could give rise to a known edge. For nonlinear image irradiance equations 
an edge constrains the possible surfaces to a small number. However, if there exists a 
surface which contains a vertex and edges emanating from it, it is unique. 

Then we discussed how singular points of an eikonal equation constrain its possible 
solutions. In particular we proved that there exists a unique (up to translation in 
the 2-direction) positive convex surface which satisfies an eikonal equation in some 
neighborhood of a singular point. 

Finally we investigated images in which b-silhouettes can be detected. An image 
irradiance equation which describes the relationship between the gradient of a surface 
whose image contains a b-silhouettc, and the shading of this surface, can be singular. 
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We showed that singular image irradiance equations may be wived using the same 
methods applied to find the integral surfaces of a continuous image irradiance equation. 
However, our ultimate ambition was to answer the following question: 

Is there a set of constraints which assure that if an image irradiance equation 
has a solution, it is unique? 

We answered this question affirmatively in chapter V. It was shown there that if three 
constraints are known to hold, the information about the imaging situation and the 
surface as captured by an image irradiance equation, allow one to reconstruct the shape 
of the surface in a unique manner. Furthermore one can easily check whether an image 
irradiance equation satisfies these constraints. It is surprising that our uniqueness 
theorem holds only when the b-silhouette is a closed curve (constraint C2). 

In order to evaluate the usefulness of our uniqueness theorem, we need to know 
which of the commonly arising image irradiance equations actually obey the above 
mentioned restrictions. In his paper on hill-shading [HOT9], Horn discusses eighteen 
different reflectance maps which are applied to solve that problem. Constraint CI holds 
for five of those reflectance maps. These equations are of the form: 

i?(p,<7) = /(p 2 +9 2 ) Cvn.i.i) 

where / is a bijection. A reflectance map of the form (VII.1.1) describes, for instance, 
the situation when the object is Lambertian and the light source and the viewer have 
the same position. Furthermore eikonal equations can be also used to automatically 
analyze images taken by a scanning electron microscope which seems to be one of the 
prime applications of our uniqueness theorem. 

There are several issues not treated here which would be interesting to investigate 
further: mutual illumination, shadows and specularities, for example. Another open 
problem is to determine the class of functions R = R{p, q) which are reflectance maps. 
Our only restriction upon these functions has been that they be C\. If a smaller set 
can be found which is known to contain all reflectance maps, better methods could 
be developed to solve the shape from shading problem in practice. A further research 
problem is to find a good numerical shape from shading algorithm. 
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Appendix 



Mathematical Details 



This appendix contains a summary of results known about first order partial 
differential equations. A thorough treatment of this material can be found in ICOH162D, 
pp.62-153] and [SMI68b, pp.259-301]. 

A.1. Basics 

For simplicity of exposition, we will discuss only first order partial differential 
equations involving a function z of two variables, x and y. The results can be generalized 
to functions of n variables in a straightforward fashion. Let p and q denote the first 
order partial derivatives of z with respect to x and y respectively. Then the relation: 

F(x,y,z,p,q) = (A- 1 - 1 ) 

where F is a function of x, y, z, p and q, is called a first order partial d ; fferent ^51 at ^ 
(abbreviated in the following by FOPDE). The relation (A.l.l) is a linear FOPDE if it 
is linear in p and q with coefficients depending only on x and y and (A.l.l) is quasi-imear 
if it is linear in p and q with coefficients depending on x,y and z. 

A function z{x,y) is called a solution to (A.l.l) if in some region of the x-y plane 
the function and its derivatives identically satisfy the equation in x and y. Such a 
function is also called an integral surface. 

The general solution to a FOPDE is a whole set of solutions, each of which satisfies 
the equation. Given a FOPDE, what kind of constraints can be imposed so that there 
is only one integral surface which satisGes both the constraints and the equation? Such 
constraints are for example boundary conditions or initial values. In section A.4 we will 
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state precisely what comprises a general solution and in sections A.6 and A.7 discuss 
what kind of constraints are necessary to pin down a particular solution. 

Unless otherwise stated, we assume that F,z and all relevant derivatives are 
continuous. 

A. 2. The Quasi-Linear FOPDE 

We consider quasi-linear FOPDE's first as their geometric interpretation is rather 
clear and so the relevant method for solving them can be explained and understood 
easily. In this case we rewrite the relation (A.1.1) as: 

a{x, y, z)p + b{x, y, z)q = c{x, y, z). (A.2.1) 

Furthermore we assume that: 

a 2 + &V0. (A.2.2) 

Suppose that the solutions to (A.2.1) are written implicitly as: 

G{x,y,z) = 0. (A.2.3) 

Differentiating (denoted by subscripts) (A.2.3) with respect to x and y yields: 

G x + G z z x = and G y -f G z z y = (A.2.4) 

or equivalently: 

p = -^ and 9 = -^- (A-2.5) 

Using these equations in (A.2.1) we obtain: 

a{x, y, z)G x {x, y, z) + b{x, y, z)G y {x, y, z) + c{x, y, z)G s (x, y, z) = 0. (A.2.6) 

Note that, in general, (A.2.1) is a nonlinear FOPDE for the function z(x,y) whereas 
(A.2.6) is a linear FOPDE for G(x, y, z). We can interpret the coefficients a, b and c 
in (A.2.6) as the components of a vector field which is defined by £ = £(x, y, z) = 
[a(x,y,z),b(x,y,z),c(x,y, z)]. Then we can rewrite (A.2.6) as: 

(£,VG) = Q (A.2.7) 

where VG denotes the gradient of G and (, ) the inner product of two vectors. Equation 
(A.2.7) expresses the fact that £ is perpendicular to VG. Since the vector VG is 
perpendicular to the surface defined by G(x, y,z) = at each point (x, y, z), we deduce 
that £ lies in the tangent plane of this integral surface at that point. 
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A field line is a curve whose tangent at every point has the same direction as the 
field vector there. An integral surface of (A.2.7) can be built up from field lines (called 
characteristic curves in this context) of the vector field f . To reiterate the previous 
statements: The tangent at each point of a characteristic curve has the same direction 
there as the vector £ and therefore, by virtue of (A.2.7), is perpendicular to the surface 
normal to the integral surface G(x, y, z) = 0. This does not mean that each quasi-linear 
FOPDE has a single solution. Such a FOPDE only constrains the possible orientations 
of the tangent planes at each point to a one-parameter manifold. As (A.2.1) is linear 
in p and q at every point of any integral surface, all feasible tangent planes intersect in 
a line which is called the Monge axis. 

We describe now a method for finding characteristic curves. Such curves can be 
defined as functions of one parameter s: x = x{s), y = y(s) and z = z(s). The vector 
[x{s),y(s),z(s)} is denoted by X - Then & = [^, ^, *£*] has the same direction 
as £ and so the outer product of ^ and £ must equal zero: 

6^-^ = 
ds as 

c ^_ o ^ = (A.2.8) 

as as 

dy dx 

a- o— = 0. 

ds ds 

The relations (A.2.8) are normally written as: 

dx:dy:dz = a:b:c. (A.2.9) 

The solutions to the equations (A.2.9) comprise a two-parameter family of curves 
in space (a family of characteristic curves), yet only a one-parameter subset of them 
generates the solutions to the FOPDE. To find this subset, an arbitrary function 
between the two free parameters obtained when solving (A.2.9) is introduced. Hence the 
general solution to (A.2.1) contains this arbitrary function. So each surface produced 
by a one-parameter family of characteristic curves is an integral surface. Conversely 
we now show that each integral surface is generated in such a fashion. 

On each integral surface z = z(x, y) the equations: 

^=a(*,y,z) ^=&(*,y,*) (A.2.10) 

ds ds 

define a one-parameter family of curves: x = x(s),y = y{s),z = z(x(s),y{s)). Note 
that on any curve in this family: 

^ = c(x,y,z) (A.2.11) 

ds 
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dz dzdx dzdy ,. -^ 

— = — — + — — = ap + bq = c. (A.2.12) 

ds ox ds ay ds 

Thus every integral surface is swept out by a one-parameter family of characteristic 
curves. The following example illustrates these ideas. 



Example: 

The following FOPDE is to be solved: 

F(x,y,z,p,q) = xp + yq — z = 0. (A.2.13) 

As a[x, y, z) = x, b(x, y,z) = y and c(x, y, z) = 1, the equations for the characteristic 
curves (A.2.9) are: 

dx:dy:dz = x:y:l (A.2.14) 

and have as their solution the two-parameter curves in space: 

y = Cix (A.2.15) 

z = C%x 

where C\ and Ci are constants. Any integral surface can be built from the curves 

described by equations (A.2.15) and every such surface is a one-parameter manifold of 

these curves. Each of these manifolds is determined through coupling C\ and Ci by 

an arbitrary relation w: 

- = C 2 = wlCx) = w(-). (A.2.16) 

x x 

Thus the solution to the FOPDE is: 

z = w{¥-)z. (A.2.17) 

x 

Writing this in parametric form gives: 

y = CiX and z = w{d)x. (A.2.18) 

A3. The General FOPDE 

We can apply methods similar to those developed in the previous section to solve 
a general FOPDE: 

F(x,y,z,p,q) = (A.3.1) 
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where we require that: F 2 , jr* ^ . ( A - 3 ' 2 > 

Our (oal is to transform the problem of Ending a solution to (A.3.1) into the problem 
Sa«»g -" of ordinal differential equations. Again, geometrical reading wd. 

^XKTjX- an integral surface. Then the quan«ie » p and * which 

P /r Ex dTA3. s ^quatZ for p and ,. To write this equation iu par^etac 
f„rouly one parameter is'needed.) The envelope of the taugent planes is a omca^ 
surface which can have several sheets and is called the Monge cone (A conical surface 
sprXed by moving a stra,ght hue which is fixed at one point, along a curve.) The 
r C !r„n» Wre refer merely to a suitable small range of tangent planes, e.g., a 
7™£TlTon£c™:tLe , can be expressed as a single-valued different.ab e 
portion ol a. sbe et M * genera tor of this cone represents a possible 

Srn of ?he' C ta n™t p'la^fit P and t caned a charade — . Jta*. 
integral surface has to fit into the field of Monge cones, l.e., has to always be tange 



Recall now that in the quasi-linear case, a Monge cone degenerates to a Monge 
axis To determ ne the solution in that case, we find the characteristic curves which 
" vel potrnave as their tangent direction the direction o the "»«»«-*£ 
Antral surface is swept out by these curves. In the case of a general FOPDE, U>e 
^meUs idea w: rk s. Ag'ain, weLve to find ^u-es which at every po^h- 

th eir tangent direction a ^^T-^R^lZXTZ^^ TZl 

words, the functions *MWfW«l^^%SZZl P 
curves determine i(s), y(s) and z(s). Yet, ^a.j.ij gives, uujr This eauation 

and a and so in order to determine p and q, another equation is needed. This equation 
T»n he obtained bv requiring that focal curves be embedded in an integral surface. (A 
^J^w2H in some neighborhood of the projection of this curve onto 
Z r « olane , is a single-valued, twice continuous differentiable function of x and 
^^t^Vh^Si-y this last condition are called ^^^ "* * 
one-oarameter family of characteristic curves sweeps out an integral surface. 

The prob em which actually has to be solved is that of finding an int egral surface 
So the soktion to this proceeds in the opposite direction from that described m the 
p eceling ^^parTgraphs. First, a set of equations, called the charactensfcc equa^on has 
to be found. The characteristic curves are a subset of the solutions of this set from 
thfchn integral surface can be built up. In the following paragraphs the technical 
prerequisites needed to find an integral surface of (A.3.1) are developed. 
P As a first task we determine the equation of the Monge cone ^F 
fixed point which has the coordinates (x,y,z). Then p and q , which satisfy (A.3.1), 
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are.written as functions of a parameter u and all feasible tangent planes at (x, y, z) are 
expressed as: 

(Z - z) = {X - x)p(u) + (r - y)q{u). (A.3.3) 

The envelope of the planes described by the previous equation defines the Monge cone, 
a conical surface whose vertex is (x,y, z). The equation of the Monge cone is found by 
eliminating u from (A.3.3) and from (A.3.4) which is obtained by differentiating (A.3.3) 

with respect to u: 

Q = {X-x)dp {Y-y)dg (A 3 4) 

du du 

Differentiating (A.3.1) with respect to the parameter u gives: 

= ^+fA (A3.5) 

v du du 

Assuming that neither ^ and gj nor F p and F q vanish simultaneously, we derive the 
following equation from (A.3.3), (A.3.4) and (A.3.5): 

^Zl = ^U. = Z ~ z . (A.3.6) 

F p F q P F p + qF q 

By substituting all possible values for p and q (i.e., all values for p and q which 
satisfy (A.3.1)) we obtain all generators of the Monge cone at the point (x, y, z). "Space 
curves having a characteristic direction at each point shall be called focal curves" 
[COHI62b, p.76]. Therefore focal curves have to satisfy the following differential 
equations: 

£ = F p 4«F f % = pF p + qF q . (A.3.7) 

da v * ds * ds 

To show that every integral surface z can be generated in such a fashion, let z = z(x, y) 
be an integral surface on which p and q is also known. Then the equations: 

d ± = Fp ^ = F q (A.3.8) 

d* da 

define a one-parameter family of curves. On these curves: 

l^p^+g* (A.3.9) 

ds ds ds 

holds and using (A.3.8) in (A.3.9) we obtain the third equation of (A.3.7): 

*E= pFp + q F q . (A.3.10) 

ds 
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The above equation is called the strip condition. "It states that the functions x{s), y{s), 
' z(s) p(s) and q{s) not only define a space curve, but simultaneously a plane tangent to 
it at every point. A configuration consisting of a curve and a family of tangent planes 
to this curve is called a strip" [COHI62b, p.77]. 

Equation (A.3.9) states that the curves defined by (A.3.7) are focal curves. Now it 
is also required that a focal curve be embedded in an integral surface. Differentiating 
(A.3.1) with respect to x and y we obtain: 

F pPx + F q q x +F z p + F x =0 (A.3.H) 

F P p y + F q q y + F a q + F v = 0. 

Using (A.3.7) and the fact that p y = q x leads to the following equations: 

i -*£+*,* -**; + **■, (A.3.12) 

ds ds da 

d Q dx i d y « rr _i_„ JP 

-T = 9x^r + 9y^" = Py F P + Wi- 
de ds ds 

Finally, using (A.3.11) in (A.3.12) yields: 

* + F.p + F. = ( A - 3 - 13 ) 

ds 

%- + F z q + F y =0. 
ds 

In summary: If a focal curve is embedded in an integral surface then along this 
curve the coordinates x, y, z and the quantities p and q satisfy the following five ordinary 
differential equations: 

d ± = F v %=F q ^= P F p + qF q (A.3.14) 

ds ds ds 

* = -(pF. + F„) f s =-(qF, + F y ). 

We can consider this system of differential equations (which are called characteristic 
equations) by itself, i.e., disregarding that we obtained it with a given integral surface 
in mind. Note first that F{x,y,z,p,q) is constant along each solution to the system 
(A.3.14) since: 

i; =F 'Ts +Fq T S ^ F 'ds^ Fx ds^" y ds 

= - F P ( P F 2 + F x ) - F 9 (qF. + F y ) + F z {pF p + qF q ) + F X F P + F y F q (A.3.15) 
=0. 
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Thus F{x,y,z,p,q) = c, where c is a constant, is a solution. to (A.3.14). The sys- 
tem of characteristic equations defines a four-parameter family of curves. By impos- 
ing the additional constraint that the solutions to (A.3.14) also satisfy the FOPDE 
F{x,y,z,p,q) = we obtain a three-parameter subset of the solutions. "Every solu- 
tion of the characteristic differential equations which also satisfies the equation F = 
will be called a characteristic strip; a space curve x(s), y{s), z{s) bearing such a strip is 
called a characteristic curve" [COHI62b, p.79]. The fact that a one-parameter subset 
of the three-parameter family of curves sweeps out the integral surface has already 
been established. Thus the problem of finding a solution to a FOPDE is equivalent to 
integrating the system of five (or equivalently four if the equations are not written in 
parametric form) ordinary differential equations (A.3.14). Note that since the charac- 
teristic curves depend on the solution, their range of influence cannot be determined in 

advance. 

In the next section we discuss the notion of a complete integral and then show how 
to choose the appropriate one-parameter subset of the solutions to the characteristic 
equations. 

A.4. General Solution, Complete and Singular Integral 

We showed in the previous section that each solution to a general FOPDE is swept 
out by a one-parameter family of curves. Thus we can write the equation of an integral 
surface as a function of the coordinates x and y and an arbitrary function of one 
variable. Such an equation is called the general solution to a FOPDE. 

Suppose now for a moment that a solution z = <P{x, y, a, b) to a FOPDE is known, 
where 4> depends on the two parameters o and 6. Then &(x, y, a, b) is called a complete 
integral if A, defined by: 

A = <Px«#yb — $xb#ya ( AA1 ) 

is not equal to zero. This condition assures that # really depends on a and 6, i.e., that 
there is no 7 = g(a, b) such that #(x, y, a, b) = $(x, y, 7). 

From the two-parameter family of surfaces defined by $(x,y,a,b), we can chose 
a one-parameter subset by introducing an arbitrary function w which relates a and 
b, e.g., by setting b = w{a). Note that the family #(x, y, a, w(a)) is a solution to the 
FOPDE. The envelope of the family <P(x, y, a, w{a)) is again a solution to the FOPDE 
since at each point it touches a member of the family ${x, y, a, w{a)), i.e., a solution. 
We obtain the equation of this envelope by eliminating the parameter a from the two 
equations: 

z = #(x,y,a,w(a)) (A.4.2) 

$ a (x, y, a, w{a)) + # 6 (x, y, a, w{a))w'(a) = 0. 

We assume throughout that all eliminations are possible and that during the course 
of this process only functions with continuous derivatives are obtained. Eliminating a 
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from (A.4.1) yields an expression involving an arbitrary function w of one variable 
which is a solution to the FOPDE and therefore the general solution. We demonstrate 
this fact analytically by differentiating the first equation of (A.4.2) with respect to x 
and y: 

2x = #x + (#a+<*W(«)K ( AA3 ) 

Zy = <p y + (# a + 4> b w'{a))a y . 
Recall that <P{x, y,a, w{a)) is a solution to the FOPDE for any choice of the parameter 
a. Using (A.4.2) (i.e., <P a + <P b w'{a) = 0) in the previous equations establishes the fact 
that for all x and y the values of z, z x and z y are the same as those of <?, $ x and # y . 

So if a complete integral of a given FOPDE is known, we can obtain the general 
solution by differentiation and by elimination of parameters. (This latter process can 
in practice be tedious or impossible but is often not necessary since every solution to 
the FOPDE is obtained by substituting all possible values for a.) In the next section 
we will show that with the help of the characteristic equations, a complete integral can 

found. 

The general solution does not comprise all solutions to a FOPDE. The envelope 
of a complete integral, the so called singular integral, is a solution which cannot be 
obtained from the general solution. The equation of the singular integral, which does 
not contain any arbitrary elements, is found by eliminating the parameters o and b 
from the equations: 

$(x, y, a,b) = z 
*a{x,y,a,*>) = (A.4.4) 

^b(^y > o»&) = °- 

In fact, we do not have to know a complete integral of a FOPDE in order to find 
the singular solution. Note that for a complete integral <P, F{x,y,4>,<P x ,<P y ) vanishes 
identically for all choices of the parameters a and b. Differentiating the FOPDE with 
respect to a and b we obtain: 

F*<P a + Fp* X a + F q $ ya = (A-4.5) 

F4,<P b + F p <P xb + F q <P yb =0. 
As <? is a complete integral, A = <P xa <P y b - $ x b$ ya is not equal to zero. Furthermore 
<P a and <P b are zero (equations (A.4.4)) and therefore we may derive the equation of the 
singular integral by eliminating p and q from: 

F p = F q =0 F = 0. (A.4.6) 

Note that equations (A.4.6) are derived from (A.4.5). (Remark: In this case we do not 
assume that F 2 + F 2 q ^ as we did to obtain the characteristic equations.) 

If a FOPDE does not contain the function z{x, y) explicitly, then no singular 
solution can exist as in this case the complete integral is of the form [COHl62b, p.24]: 

z = <P{x,y t a) + b (A.4.7) 

and the condition <P b — cannot be fulfilled. 
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A.,5. A Method for Finding the Complete Integral 

In the previous sections we showed two methods for finding the solutions to a given 
FOPDE. Here, we explain how to find the complete integral using the characteristic 
equations. Furthermore we describe how to determine a one-parameter subset from the 
four-parameter family of solutions to the characteristic equations. 

First let us discuss a special form of a FOPDE called PfafTs equation: 

/(x, y, z)dx + g{x, y, z)dy + h{x, y, z)dz = 0. (A.5.1) 

In the case where h = and / and g depend only on x and y, this equation degenerates 
to an ordinary differential equation called an exact differential equation: 

f(x, y)dx + g{x, y)dy = 0. (A.5.2) 

The equation is called total if / and g satisfy the integrability condition: 

f y {x,y) = g x (x,y). (A5.3) 

In the case of a total differential equation it is easy to find a solution to (A.5.2). On 
each simply connected region we can determine a function H (x, y) such that: 



dH „, , , dH 

Then: 



^ = /(x,y) and ^- = ff(a:,y). (A.5.4) 

dx oy 



dH = /(x, y)dx + g{x, y)dy (A.5.5) 

and the equation dH = are both equivalent to (A.5.2). Thus H (x, y) = c, where c 
is a constant, is a solution to (A.5.2) and the function H can be found by integrating 

(A.5.4). 

In the case where a FOPDE is exact but not total, an integration factor /x(x, y) 
can be always introduced such that the equation: 

(ifdx + l^gdy = (A.5.6) 

is total, i.e., such that (/i/) v = {jig) x . Equivalently, /x(x,y) has to be a solution to the 
following FOPDE which we can solve using the method of characteristic curves: 

M(/y-<fc) + Mv/-M*0 = O. (A.5.7) 

Equation (A.5.1) is also easy to solve if its left hand side is a total differential of a 
function H (x, y, z) i.e., if: 

. dH dH dH , A . 

/== ax" g= ~dy- k =-dz- (A - 58) 
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. For the previous equations to hold, i.e., for the functions ■/, g and h to be the first 
order partial derivatives of a function H, it is necessary that the following integrability 
conditions be satisfied: 

dy dx dz dy dx~dz' ^'^ 

In a simply connected region these integrability conditions assure the existence of a 
function H(x, y, z) which satisfies (A.5.8) and can be calculated as: 

H{x, y, z) = / (fdx + gdy + hdz) + C (A.5.10) 

« / (xo,yo,*o) 

where (x ,y Q ,z ) is a fixed point and C is a constant. Clearly H(x,y,z) = c, where c 
is a constant, is a solution to (A.5.1). 

In (A.5.9) is not satisfied, it is again desirable to find an integration factor /i(x, y, z) 
such that the expression nfdx + ngdy -f fxhdz is a total differential of a function. In 
contrast to the case of Pfaff's equation in two variables, it is not always possible to find 
such a factor. For such a // to exist, it is necessary that the following equation holds: 

f{9, ~ h y ) + g(h x - f M ) + h{f y - g x ) = 0. (A.5.11) 

It can also be shown (simple, but tedious) that in a simply connected region the 
previous equation is sufficient for (A.5.1) to possess a one-parameter family of solutions 
H(x, y, z) = c. We demonstrate now how to construct such a function H(x, y, z). 
First consider the abbreviated equation: 

f{x, y, z)dx -f g(x, y, z)dy = 0. (A.5.12) 

This is a Pfaff's equation in the two variables x and y, with z as a parameter. Thus a 
solution to it can be always found (if necessary with the help of an integrating factor 

$(x, y) = u{x, y,z) = c (A.5.13) 

where c is a constant. Note that: 

x , du , du 

Xf== di and Xi7== ay"- (A>514) 

Now we define a function S which depends upon the three variables x, y and z: 

S(x, y,z) = \h——. (A.5.15) 
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Using (A.5.13) we redefine 5 as a function of x,u and z, say: 

T(x,u,z) = S{x,y,z). (A.5.16) 

Suppose that T is independent of x. Then we can find H(x,y,z) by solving another 
Pfaff's equation in the variables u and z. In fact, it is easy to see that T is independent 
of x. One just has to prove that |£ = which we do using equations (A.5.14) and 
(A.5.15): 

±(X ft -S) = u„ = ±(X/) 

^(\h - S) = u y , = £(Xg) (A.5.17) 

£(X 5 ) = „, s = A (X/) . 

Equations (A.5.17) written out in full are: 

S x = h\ x - f\ z + \{h x - f z ) (A.5.18) 

-S y = gX a - h\ y + X((7 a - h y ) (A5.19) 

= /X y - 9 X x + X(/ 1/ -g I ). (A.5.20) 

Multiplying (A.5.18) by g, (A.5.19) by /, and (A.5.20) by h and then adding up the 
three equations using condition (A.5.11) yields: 

gS x — fS y = 0. (A.5.21) 

Differentiating (A.5.16) with respect to x and y gives: 

Sx = T x + T u u x (A.5.22) 

By combining (A.5.14), (A.5.21) and (A.5.22), we obtain: 

= gS x - fS y = gT x + gT u u x - fT u u y = (A.5.23) 

= 9T X + \fgT u - \fgT u = gT x . 

Because g^Owe may conclude that T x = 0. Thus we rewrite (A.5.16) as: 

S{x,y,z) = T(u,z). (A.5.24) 

Equation (A.5.1), after multiplication by X and using the expressions (A.5.14) and 
(A.5.15) becomes: 

\{fdx -}- gdy + hdz) = u x dx + u y dy + (u z + T)dz = (A.5.25) 



97 

or, equivalently: 

du + T{u,z)dz = 0. (A.5.26) 

This is again a Pfaffs equation in two variables, which can always be solved using an 
integrating factor. Its solution is ip(u, z) = c where c is a constant. Therefore the 
solution to (A.5.1) is: 

H{x, y, z) = ip{u{x, y, z), z) = c (A.5.27) 

which is a one-parameter manifold. 

Finally, we describe a method for finding a complete integral of a general FOPDE 
F(x,y,z,p,q) = 0. The basic idea is that the total differential of a solution z(x,y) to 
the FOPDE: 

dz = pdx -f- qdy (A. 5. 28) 

can be interpreted as Pfaff's equation in the variables x, y and z. In order to do so, p 
and q must be expressed as functions of x, y and z. So assume that two functions / 
and g can be found such that: 

P = f{x,y,z,a) and q = g(x, y, z, a) (A.5.29) 

where a is an arbitrary constant. Then the FOPDE and (A.5.11) (with h = — 1): 

f9z-9f,-fy + 9 x =0 (A.5.30) 

are satisfied. In this case the solution to (A.5.12) is a one-parameter manifold, but since 
a parameter o is already build into the equation, this solution contains two parameters 
and is the complete integral. Hence the problem of finding the appropriate functions 
/ and g has to be solved. 

Suppose that a function G(x, y, z, p, q) exists such that the following equation can 
be solved for p and q (or equivalently for / and g): 

F{x,y,z,p,q) = (A.5.31) 

G(x,y,z,p,q) = a. 

For this to be possible, the inequality F p G q — F q G p 7^ has to hold. We have to 
prove that / and g obtained in such a fashion satisfy (A.5.30) identically in the three 
variables x, y and z. Differentiating the equations (A.5.31) with respect to x, y and z 
gives: 

F x + p x F p + q x F q = G x + Px G p + q x G q =0 

F v +PyF p -\-q y F q = G y + Py G p + q y G q = (A.5.32) 

F z + p z F p + q z F q = G z + p z G p + q z G q = 0. 

After expressing p z , q z , p y and q x from the previous equations in terms of the derivatives 
of F and G and substituting these expressions into (A.5.30) we obtain a linear FOPDE 
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for the function G: 

F P G X + F q G y + (pF p .+ qF q )G, - (F x + P F,)G P - {F y + qF,)G q = 0. (A.5.33) 

We can solve this equation by the method of characteristic curves. The appropriate 
system of characteristic equations is the same as the one for F{x,y,z,p,q) = 0: 

dx:dy:dz:dp:dq - F p :F q :{pF p + qF q ): - [pF, + F x ): - {qF, + F v ). (A.5.34) 

Only one solution to these equations which is independent of F and contains at least 
one of the variables p and q is needed. Such an integral is the desired function G 
and will always exist since the solution to the characteristic equations comprises a 
four-parameter family: 

v i {x,y,z,p,q) = C i t = l,2,3,4. (A.5.35) 

The V{ are independent and at least one of them must contain either p or q. 

The method we have described is due to Lagrange and Charpit. It has the 
advantage over the method of characteristic curves discussed in section A.4 in that 
one needs only to find a single integral of (A.5.34) rather then a four- parameter family 
of curves. 

A.6. The Initial Value Problem for Quasi-Linear FOPDE's 

In this section we attack the problem of how to determine a particular solution of 
a FOPDE, once the general solution is known. First we consider a quasi-linear FOPDE: 

a{x,y,z)p + h{x,y,z)q = c{x,y t z). (A.6.1) 

In particular we discuss the problem of how to find an integral surface z{x, y) of (A.6.1) 
which contains a given curve C in space. In the literature, e.g., [COHI62b, p. 40], this 
problem is referred to as Cauchy's problem. Clearly the following questions have to be 
answered: 

1) What conditions on C are necessary to make this problem is solvable? 

2) When is such a solution unique? 

Let C be given by three continuous differentiable functions of a parameter t: x{t), y{t), z(t). 
Furthermore, we assume that the projection of C onto the x-y plane (referred to as C ) 
does not contain double points and that x\ + y\ ^ 0. If C contains double points, an 
integral surface with self-intersections is obtained, hence, z is not everywhere a single- 
valued function of x and y which implies that along the line of intersection, p and q 
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are discontinuous. Now to find a solution to the FOPDE containing C, we construct 
a characteristic curve through each point of C. The equations for the characteristic 
curves depend upon two parameters: 

x = x(s, t) y = y(s, t) z = z(s, t). (A.6.2) 

Note that the functions x, y, z are still continuously differentiable. To get the equation 
of the integral surface, we must eliminate the parameters s and t from the previous 
equations, i.e., we must express s and t in terms of x and y. A sufficient condition 
for this to be possible is that the functional determinant A, which is specified by the 
following equation, does not vanish along the curve C: 

as at as at 

Using the characteristic equations, we rewrite (A.6.3) as A = a-j£ — b^ft. Thus if 
A 7^ 0, we may express z as a function of x and y and be assured that C lies on the 
surface. This solution is also unique which is a consequence of the following lemma: 

Lemma [COHI62b, p.64]: Each characteristic curve which has one point in 
common with an integral surface lies completely on this surface. 



This lemma follows from the uniqueness theorem for solutions to ordinary differential 
equations. 

We can interpret the determinant A as the outer product of the two vectors: 

«,-(*) »d 6-(*) (A.6.<) 

which are the projections of the tangent and the characteristic direction onto the x-y 
plane, respectively. In the special case where A vanishes along C, these two directions 
coincide and we may deduce that C has one of the following three properties: 

1) C is a characteristic curve. 

2) C is the envelope of the characteristic curves (called an edge of regression). 

3) Co is the envelope of the projections of the characteristic curves onto the 
x-y plane. 

We discuss when case 1 occurs first. From A = we infer that: 

~ = I*- (A-6.5) 

adt bdt K J 
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Using x(t) and y(t) (from the equation for C) in z(x, y) the following equation has to 
hold along C: 

dz dx , dy 

Tt =z *Tt JrZ Ut =ap + bq==c - < A - 6 - 6 ) 

This means that C satisfies the characteristic equations and is therefore a characteristic 
curve. In this case, there exist infinitely many surfaces through C which satisfy the 
FOPDE. To see this, just choose another curve C which has a point P in common with 
C. Now to construct the solutions through C, a characteristic curve is passed through 
every point of C, in particular through P. The characteristic curve through P is C, thus 
an integral surface through C contains C. In this manner we can construct infinitely 
many integral surfaces which contain C. They all meet along the characteristic curve 
C which therefore can be viewed as a branch curve. 

One assumption made throughout should be stressed again here: we require the 
solutions to a FOPDE to be continuous and continuously differentiable in some neigh- 
borhood of C. It might be possible to find a solution z through C, along which A 
vanishes, without C being a characteristic curve as occurs in cases 2 or 3 mentioned 
above. However, the derivatives of z are then not continuous along C. This fact is 
illustrated in the following example. 

Example: 

We wish to solve the following equation: 

F{x, y, z, p, q) = 3(z — yfp — q = 0. (A.6.7) 

The characteristic equations are: 

dx 

ds 
dy 



Ts = 3 <'-'> 2 



dz 
ds 

The solution to these equations, with the initial values XQ,y 0i z Q , is: 

x = (z — y + s) 3 + x Q — (z — j/o) 3 

y = —s + y (A.6.9) 

z — z . 

Now we impose the constraint that the integral surface passes through C which is given 
by: 

i = y = t z = t. (A.6.10) 
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Note that C is not a characteristic curve. Setting the initial values x , y , z (i.e., x, y, z 
for s = 0) to: 

x = y = t z = t, (A.6.11) 

(A.6.9) becomes: 

x = s 3 

y = —s + t (A.6.12) 

z = t. 

In this case the determinant A is: 

dxdy dxdy 2 

ds dt dt ds ' \ ■ • ) 

Thus along the curve C (i.e., s = 0) A = 0. 

There is, however, a solution to the FOPDE which contains C: 

z = i* + y. (A.6.14) 

Note that p = ^x~3 does not exist along C (as x = there). 

In the case of a linear FOPDE we can make some further statements about the 
solution in the case where A vanishes along C. Here, the integral surfaces are cylindrical 
surfaces perpendicular to the x-y plane, i.e., the function defining an integral surface 
is independent of z. The linear FOPDE is: 

a{x, y)p + b(x, y)q = c(x, y). (A.6.15) 

Recall that A is defined as: 

A = x a y t — x t y a . (A.6.16) 

Using the characteristic equations: 

x a = a (A.6.17) 

Vs = b 

the following equation is obtained for A s : 

A 5 = a s y t + ay st — b s x t — bx at = 

= a sVt + ab t — b s x t — a t b. (A.6.18) 

Note that if A 5 (the first order partial derivative of A with respect to s) and A vanish 
along C, then A vanishes everywhere. The proof of this last assertion follows from the 
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existence and uniqueness theorem for ordinary differential equations. By differentiating 
a and 6 with respect to s and t and using relations (A.6.17), equations (A.6.19) are 
obtained: 

a a = a x a -f a y b 

at = a x x t + a y y t (A.6.19) 

b a = 6 x o -j- b y b 

b t = b x x t -f b y y t 

and so: 

A 5 = (a, + b y )A. (A.6.20) 

To express x as a function of y and 2, say, x = f{y,z), we have to assume that 
y s zt — z s y t =£0 along C. 

Now, in order to prove that the integral surface z is a cylindrical surface, it suffices 
to show that f z = 0. Differentiating x with respect to s and t we derive the following 
equations: 

X» = fyVs + fzZ a 

*t = fyyt + f z z t . (A.6.21) 



Then: 

from which it follows that 



A = f 2 (z s y t — z t y a ) (A.6.22) 

fz = 0. (A.6.23) 

A. 7. The Initial Value Problem for General FOPDE's 

Here we pose the question: What are the constraints necessary to determine a 
solution to a general FOPDE uniquely? Clearly more information than in the quasi- 
linear case is needed as now the solutions to the characteristic equations form a three- 
parameter family of curves. So let C be a curve given by x(t), y{t) and z[t) such that 
neither C nor its projection onto the x-y plane have double points. Furthermore p(t) 
and q(t) along C have to be specified such that the condition: 

dz dx dy 

Tt= p Tt+ q Tt ( AJ - 1 ) 

holds and the FOPDE (A.3.1) (i.e., F = 0) is identically satisfied in t. Thus the 
functions x(t),y(t),z(t),p(t) and q(t) define an initial strip which we denote by C x . 
From now on, the procedure is very similar to the one used in solving the initial value 
problem for a quasi- linear FOPDE. Through every element of d a characteristic strip 
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is constructed which can be written as x(s, t), y(s, t), z(s, t),p(s, t) and q(s, t). To be able 
to express the parameters s and t in terms of x and y: 

A = x a y t — x t y 3 = F p y t — F q x t (A.7.2) 

cannot vanish identically along the initial strip. Assuming this holds, z,p and q can be 
expressed in terms of x and y. It then remains to check whether p and q written in 
such a fashion are the first order partial derivatives of the integral surface z(x, y). This 
involves showing that the quantities U and V: 

U = zt—pxt— qy t (A.7.3) 

V = z a — px 3 — qy a 

vanish identically. Since we assumed that A^Owe deduce from the previous equations 
and: 

= z t — z x x t — z y y t 

= z s — z x x a — z y y a (A.7.4) 

that z x = p and z y = q. Recall now the characteristic equations: 



ds p ds q ds 



Fv ? = *■« ™= P F P + qF q . (A.7.5) 



Using the first two of these equations in the last one, we obtain: 

— = — 4- ^L ( A 7 6) 

ds ds ds 

which implies that V vanishes identically. Now to prove that also U vanishes identically, 
note that: 

-5- — z*t — PsXt — px at — q s y t — qy at (A.7.7) 

OS 

dV 

-5- = *«t — Ptx a — px at — q t y a — qy at . (A.7.8) 

ds 

Subtracting (A.7.7) from (A.7.8) yields: 

dU dV . .. __. 

-z 37 = —PsXt — Ptx a 4- q s y t — q t y a . (A.7.9) 

OS ot 

Taking into account the characteristic equations and the fact that V = implies 
4j£ = 0, we rewrite the previous equation as: 

flTT 

-z- = PtF p + q t F q + x t F x + y t F y + {px t + qy t )F M . (A.7.10) 

OS 
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However, the FOPDE F = holds identically in a and t. Differentiating F with respect 
to t gives: 

F x z t +F v y t + F M z t +F p p t +F q q t =0. (A.7.11) 

Using the previous equation in (A.7.10) we obtain: 

*g-=^F,U. (A.7.12) 

08 

For any fixed t, (A.7.12) is an ordinary differential equation for U as a function of 8 
with the solution: 

U{8) = U{0)eJo ~ F ' U . (A.7.13) 

Since by assumption, U(0) is zero, U vanishes everywhere. 

To summarize the previous results: given a curve x(t), y(t) and z(t) along which 
p(f) and q(t) are known such that: 

dz dx dy 

~dl~ P ~dt +9 'dt 

F{x{t), y{t), z{t), p(t), q(t)) = (A.7.14) 

A = F p y t -F q x t ^0, 

there exists a unique integral surface through the initial strip. We obtain a unique 
surface because the solution to the characteristic equations is uniquely determined by 
their initial values. 

The exceptional case where A = along C\ is analogous to the one discussed in the 
previous section: there are infinitely many integral surfaces through C\ if and only if it 
is a characteristic strip. Again we can view any characteristic curve as a branch element 
since, on either side of it, there can be another member of the family of solutions to a 
FOPDE, while along such a curve the first derivatives are continuous. Note that higher 
order derivatives along the curve may be discontinuous. If C% is only a focal strip along 
which A = 0, then it might be possible to find an integral surface z containing it. As 
in the quasi-linear case, this surface does not have continuous derivatives. 

Finally, suppose C degenerates to a point P with coordinates (xo, yo, zq). Then 
the strip condition is identically satisfied for all po and qo which satisfy the FOPDE, 
i.e., for all po and go which determine the feasible tangent planes at P. So po and qo 
can be written as functions of a parameter t. If the quantities Xo,yo,z ,po(t) and qo{t) 
are used as initial values when solving the characteristic equations, a unique integral 
surface which has a conical singularity at P is obtained. It is called the integral conoid 
of the partial differential equation at P. 
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Appendix II 



Transformation 



In this appendix we show that any solution to an image irradiance equation of the 
form (B.l) can be obtained from the solutions to an image irradiance equation of the 
form (B.5). Let 

f{Ap 2 + 2Bpq + Cq 2 + 2Dp + 2Eq)=E{x,y) (B.l) 

be an image irradiance equation where / is a bijection and A, B, C, D and E are real 
constants such that 6 > and AS < 0, where 6, A and 5 are defined by: 



'6= AC — B* 
A B D 

A = B C E 
DEO 

S=A + C. 



(B.2) 



As / is a bijection, equation (B.l) and the following transformed equation hare the 
same solutions: 



V + 2BP9 + Cq 2 + 2Dp + 2Eq = f-\E{x,v)). 
Equivalently, we can write the previous equation as: 

Ap 2 + 2Bpq + Cq 2 + 2Dp + 2Eq = E{x, y). 



'(B.3) 



(B.4) 
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Now let z = z(x, y) be a solution to: 

p 2 + g 2 = E(x,y) + c (B.5) 

where x, y and c are defined by: 

B , y/AC-B 2 

x = x H y 

y = >/Cs (B.6) 

c = Aa 2 + 2Ba/3 + C0 2 

and a and /? are defined by: 

CD — BE 



a = 



(B.7) 



AC — B 2 
AE-BD 

P ~ AC-B 2 ' 

Then z(x, y) = z(x(x, y), y(*> £)) — ax — /?y is a solution to (B.3). 

Proof: The proof proceeds by diagonalization of a quadratic form. First we express x 
and y as functions of x and y: 



x =■ 



y (B.8) 



Cx — By 

v = 



^/C{AC — B 2 ) 
The first order partial derivatives of z are abbreviated by p and q: 

- as 



9 ~dy" 
In the following equations we express p and q in terms of p and g and vice versa: 

V = Q—, = — <* 
\MC — B 2 

\/C ^(AC - B 2 ) 



j?(p + Q ) + C(g + /?) 
9 = — (P + <*)• 



IM-i i -iimJIP 
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From this last eqa&yon «e dedttce tlui|; 

^ , i - ,- j™ ft 

which ia the same aa (B.4). | — ---- 
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